NAVAL 

POSTGRADUATE 

SCHOOL 

MONTEREY,  CALIFORNIA 


THESIS 


FOCUSING  ISAR  IMAGES  USING  FAST  ADAPTIVE  TIME- 
FREQUENCY  AND  3D  MOTION  DETECTION  ON  SIMULATED 

AND  EXPERIMENTAL 

RADAR  DATA 

by 

Wade  Brinkman 

June  2005 

Thesis  Advisor: 

Michael  A.  Morgan 

Co-Advisor: 

T.  Thayaparan 

Second  Reader: 

Jeffrey  B.  Knorr 

Approved  for  public  release,  distribution  is  unlimited 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  No.  0704-0188 
Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including 
the  time  for  reviewing  instruction,  searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and 
completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any 
other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington 
headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite 
1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project 
(0704-0188)  Washington  DC  20503.  _ _ 

I.  AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

_ June  2005 _ Master’s  Thesis _ 

4.  TITLE  AND  SUBTITLE:  Focusing  ISAR  Images  using  Fast  Adaptive  Time-  5.  FUNDING  NUMBERS 
Frequency  and  3D  Motion  Detection  on  Simulated  and  Experimental  Radar  Data 

6.  AUTHOR(S)  Wade  Brinkman  _ 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES)  8.  PERFORMING 

Naval  Postgraduate  School  ORGANIZATION  REPORT 

Monterey,  CA  93943-5000  NUMBER 

9.  SPONSORING  /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES)  10.  SPONSORING/MONITORING 
N/A  AGENCY  REPORT  NUMBER 

II.  SUPPLEMENTARY  NOTES  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official 

policy  or  position  of  the  Department  of  Defense  or  the  U.S.  Government. _ 

12a.  DISTRIBUTION  /  AVAILABILITY  STATEMENT  12b.  DISTRIBUTION  CODE 

Distribution  Statement  (mix  case  letters) _ 

13.  ABSTRACT  (maximum  200  words) 

Optimization  algorithms  were  developed  for  use  with  the  Adaptive  Joint  Time-Frequency  (AJFT)  algorithm  to  re¬ 
duce  Inverse  Synthetic  Aperture  Radar  (ISAR)  image  blurring  caused  by  higher-order  target  motion.  A  specific 
optimization  was  then  applied  to  3D  motion  detection.  Evolutionary  search  methods  based  on  the  Genetic  Algo¬ 
rithm  (GA)  and  the  Particle  Swarm  Optimization  (PSO)  algorithm  were  designed  to  rapidly  traverse  the  solution 
space  in  order  to  find  the  parameters  that  would  bring  the  ISAR  image  into  focus  in  the  cross-range.  3D  motion 
detection  was  achieved  by  using  the  AJTF  PSO  to  extract  the  phases  of  3  different  point  scatterers  in  the  target 
data  and  measuring  their  linearity  when  compared  to  an  ideal  phase  for  the  imaging  interval  under  investigation. 
The  algorithms  were  tested  against  both  simulated  and  real  ISAR  data  sets. 

14.  SUBJECT  TERMS  Inverse  Synthetic  Aperture  Radar,  Adaptive  Joint  Time-Frequency  15.  NUMBER  OF 

Algorithm,  Genetic  Algorithm,  Particle  Swarm  Optimization,  3D  Motion  Detection,  ISAR,  AJTF,  GA,  PAGES 


PSO 

143 

16.  PRICE  CODE 

17.  SECURITY 

18.  SECURITY 

19.  SECURITY 

20.  LIMITATION 

CLASSIFICATION  OF 

CLASSIFICATION  OF  THIS 

CLASSIFICATION  OF 

OF  ABSTRACT 

REPORT 

PAGE 

ABSTRACT 

Unclassified 

Unclassified 

Unclassified 

UL 

NSN  7540-01-280-5500 

Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  239-18 


1 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


11 


Approved  for  public  release,  distribution  is  unlimited 


FOCUSING  ISAR  IMAGES  USING  FAST  ADAPTIVE  TIME-FREQUENCY 
AND  3D  MOTION  DETECTION  ON  SIMULATED  AND  EXPERIMENTAL 

RADAR  DATA 


Wade  H.  Brinkman 
Captain,  Canadian  Forces 
B.S.,  Royal  Military  College  of  Canada,  1997 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  ELECTRICAL  ENGINEERING 


from  the 


NAVAL  POSTGRADUATE  SCHOOL 
June  2005 


Author:  Wade  Brinkman 


Approved  by:  Michael  A.  Morgan,  Thesis  Advisor 


Thayananthan  Thayaparan,  DRDC 
Co-Advisor 


Jeffrey  B.  Knorr 
Second  Reader 


John  P.  Powers 

Chairman,  Department  of  Electrical  and  Computer  Engineering 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


IV 


ABSTRACT 


Optimization  algorithms  were  developed  for  use  with  the  Adaptive  Joint  Time- 
Frequency  (AJFT)  algorithm  to  reduce  Inverse  Synthetic  Aperture  Radar  (ISAR)  image 
blurring  caused  by  higher-order  target  motion.  A  specific  optimization  was  then  applied 
to  3D  motion  detection.  Evolutionary  search  methods  based  on  the  Genetic  Algorithm 
(GA)  and  the  Particle  Swarm  Optimization  (PSO)  algorithm  were  designed  to  rapidly 
traverse  the  solution  space  in  order  to  find  the  parameters  that  would  bring  the  ISAR  im¬ 
age  into  focus  in  the  cross-range.  3D  motion  detection  was  achieved  by  using  the  AJTF 
PSO  to  extract  the  phases  of  3  different  point  scatterers  in  the  target  data  and  measuring 
their  linearity  when  compared  to  an  ideal  phase  for  the  imaging  interval  under  investiga¬ 
tion.  The  algorithms  were  tested  against  both  simulated  and  real  ISAR  data  sets. 
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EXECUTIVE  SUMMARY 


Inverse  synthetic  aperture  radar  (ISAR)  imaging  is  an  effective  way  to  acquire 
high  resolution  images  of  targets  of  interest  at  long  range  and  as  such  is  an  irreplaceable 
tool  in  the  task  of  non-cooperative  target  recognition  (NCTR)  of  both  ships  and  aircraft. 
The  Canadian  Air  Force  is  currently  upgrading  her  fleet  of  CP- 140  Aurora  maritime  pa¬ 
trol  aircraft  to  possess  NCTR  through  ISAR  in  order  to  increase  the  capability  of  the  Ca¬ 
nadian  Forces  in  both  sovereignty  patrols  of  Canadian  territory  and  the  protection  of  the 
Canadian  Patrol  Frigates  and  allied  ships  operating  abroad  as  a  coalition.  The  United 
States  Navy’s  equivalent  aircraft,  the  P-3  Orion,  currently  possesses  this  capability. 
Therefore,  effective  ISAR  imaging  will  have  a  real  impact  in  the  decision-making  proc¬ 
ess  in  future  military  operations  involving  both  Canadian  and  American  forces. 

One  significant  problem  with  ISAR  image  formation  is  the  assumption  of  time  in¬ 
variance  of  the  Doppler  frequency  used  to  resolve  the  image  in  the  cross-range.  Time- 
variant  Doppler  becomes  present  in  an  ISAR  signal  when  an  aircraft  is  maneuvering  or  a 
ship  is  pitching  and  rolling  during  the  coherent  processing  interval  and  is  typically  re¬ 
ferred  to  as  motion  error.  Because  of  this  motion  error  and  the  fast  Fourier  transform’s 
(FFT)  basic  assumption  of  stationary  signals,  the  use  of  the  FFT  to  resolve  the  image  in 
the  cross-range  will  cause  extensive  blurring  and  leave  the  image  unrecognizable  even  to 
the  most  experienced  ISAR  operators. 

One  proven  method  of  motion  compensation  is  the  adaptive  joint  time-frequency 
(AJTF)  algorithm.  However,  this  algorithm  is  not  without  two  significant  weaknesses  of 
its  own.  The  first  problem  is  that  the  computational  burden  of  the  exhaustive  search  used 
to  extract  the  motion  compensation  parameters  is  quite  large  which  limits  its  usefulness 
in  an  operational  situation.  The  second  problem  with  the  AJTF  algorithm  is  that  the  tar¬ 
get  must  follow  the  mathematical  model  of  the  AJTF  which,  in  particular,  states  that  rota¬ 
tional  motion  must  be  confined  to  a  2D  plane  (i.e.,  no  3D  motion). 

To  address  the  first  problem,  two  evolutionary  searches,  a  Genetic  Algorithm 
(GA)  and  a  Particle  Swarm  Optimization  (PSO),  were  investigated,  developed  and  evalu- 
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ated  in  order  to  take  advantage  of  their  properties  of  rapid  convergence.  Both  the  GA  and 
the  PSO  performed  excellently  during  their  testing  against  simulated  and  real  radar  data 
in  the  extraction  of  the  search  parameters  from  the  solution  space.  They  would  typically 
converge  to  the  search  parameters  required  for  focusing  in  less  than  50%  of  the  cost  (or 
fitness)  function  calculations  required  by  the  exhaustive  search.  Also,  the  PSO  proved  to 
be  a  much  more  efficient  algorithm  for  this  optimization  problem  with  its  performance 
consistently  exceeding  that  of  the  GA. 

To  address  the  second  problem,  the  PSO  algorithm  designed  above  was  combined 
with  the  3D  motion  detection  algorithm  developed  at  DRDC  -  Ottawa  by  Thayaparan. 
The  resulting  algorithm  was  able  to  effectively  and  efficiently  sort  good  imaging  intervals 
in  the  real  radar  data  set  from  poor  imaging  intervals  which  violated  the  2D  motion  as¬ 
sumption  of  the  mathematic  model.  It  was  shown  that  the  imaging  interval  indicated  as 
possessing  a  low  degree  of  3D  motion  created  well-focused  images  after  being  motion 
compensated.  Imaging  intervals  that  were  indicated  as  possessing  a  high  degree  of  3D 
motion  did  not  focus  even  after  being  motion  compensated.  These  imaging  intervals  can 
now  be  discarded  instead  of  presented  to  the  ISAR  operator. 

This  thesis  draws  from  many  sources  in  the  radar  signal  processing  community 
and  cites  them  appropriately.  There  are  two  contributions  made  by  this  thesis  that  are 
novel  in  their  entirety.  First,  the  automatic  detection  of  focus  using  the  FFT  for  determin¬ 
ing  basis  function  suitability  in  the  AJTF  algorithm  is  new,  effective  and  fast  and  not 
used  in  any  of  the  cited  references.  Secondly,  the  AJTF  PSO  algorithm,  with  the  PSO 
based  on  the  work  in  [7],  is  unique  to  its  application  to  ISAR  image  focusing.  The  work 
on  3D  motion  detection  draws  heavily  from  the  work  in  [12]  and  provides  two  significant 
improvements.  First,  the  measure  of  non-linearity  in  an  imaging  interval  is  based  on  the 
percentage  of  deviation  from  the  line  of  least-squares  fit  as  opposed  to  the  basing  it  on  the 
difference  between  the  line  of  least-squares  fit  and  the  actual  line  depicting  the  linearity 
of  phases.  Secondly,  the  AJTF  PSO  developed  in  the  first  part  of  this  thesis  is  applied  to 
the  3D  motion  detection  algorithm.  This  creates  an  algorithm  that  evaluates  the  degree  of 
3D  motion  in  an  imaging  interval  both  rapidly  and  accurately  which  is  useable  by  an 
ISAR  operator  in  an  operational  situation. 
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I.  INTRODUCTION 

A.  RESEARCH  REQUIREMENT  AND  OBJECTIVES 

Inverse  synthetic  aperture  radar  (ISAR)  imaging  is  an  effective  way  to  acquire 
high  resolution  images  of  targets  of  interest  at  long  range  and  as  such  is  an  irreplaceable 
tool  in  the  task  of  non-cooperative  target  recognition  (NCTR)  of  both  ships  and  aircraft. 
The  Canadian  Air  Force  is  currently  upgrading  her  fleet  of  CP- 140  Aurora  maritime  pa¬ 
trol  aircraft  to  possess  NCTR  through  ISAR  in  order  to  increase  the  capability  of  the  Ca¬ 
nadian  Forces  in  both  sovereignty  patrols  of  Canadian  territory  and  the  protection  of  the 
Canadian  Patrol  Frigates  and  allied  ships  operating  abroad  as  a  coalition.  The  United 
States  Navy’s  equivalent  aircraft,  the  P-3  Orion,  currently  possesses  this  capability. 
Therefore,  effective  ISAR  imaging  will  have  a  real  impact  in  the  decision  making  process 
in  future  military  operations  involving  both  Canadian  and  American  forces. 

As  explained  in  [1],  [2]  and  [3]  and  repeated  in  Chapter  II  of  this  thesis,  one 
drawback  of  ISAR  is  that  higher  order  target  motion,  observed  when  ships  are  pitching 
and  rolling  in  rough  seas  or  when  aircraft  are  maneuvering,  introduces  time-varying  Dop¬ 
pler  into  the  radar  signal.  It  is  this  time-varying  Doppler  that  causes  severe  image  blur¬ 
ring  in  the  cross-range  during  the  ISAR  signal  processing  stage  and  effectively  renders 
the  target’s  image  unrecognizable  even  to  the  most  experienced  ISAR  operators.  There¬ 
fore,  for  ISAR  to  realize  the  operational  goal  of  NCTR,  methods  must  be  developed  to 
compensate  for  the  higher  order  motion.  The  first  objective  of  this  thesis  was  to  continue 
the  work  by  Thayaparan  [1]  and  Chen  [3]  on  the  Adaptive  Joint  Time  -  Frequency  Algo¬ 
rithm  (AJTF)  and  to  examine  different  methods  of  efficiently  removing  time-varying 
Doppler  from  the  radar  signal  in  order  to  form  images  that  are  well  focused  in  the  cross¬ 
range. 

The  AJTF  algorithm  in  [1]  and  [3]  assumes  2D  motion  on  a  rotational  axis  during 
the  required  ISAR  dwell  interval.  Unfortunately,  true  targets  do  not  always  follow  this 
assumption  and  within  the  radar  signal  there  exists  3D  motion  which  cannot  be  compen¬ 
sated  by  the  AJTF  algorithm.  One  way  to  deal  with  the  existence  of  3D  motion  in  an 
ISAR  imaging  data  set  is  to  determine  to  what  degree  3D  motion  exists  in  an  imaging  in- 
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terval  of  the  data  set.  The  imaging  intervals  with  a  small  degree  of  3D  motion  are  imaged 
while  the  intervals  that  possess  a  high  degree  of  3D  motion  are  discarded.  The  research  in 
[3],  [4]  and  [12]  provide  the  background  for  3D  motion  detection  using  the  nonlinearity 
of  phase  method.  Therefore,  the  second  objective  of  this  thesis  is  to  apply  the  fast  AJTF 
algorithm  found  in  the  first  part  of  this  thesis  to  3D  motion  detection  methods  found  in 
[3],  [4]  and  [12], 

B.  SCOPE  AND  ORGANIZATION 

1.  Scope 

This  thesis  uses  two  ISAR  data  sets,  one  which  is  a  B727  simulated  data  set  cre¬ 
ated  by  Chen  [18]  and  is  designed  to  display  the  ideal  motion  error  found  in  a  blurred 
ISAR  image  and  is  available  from  the  Naval  Research  Laboratory  website.  An  example 
of  the  unfocused  and  focused  image  is  found  in  Figure  1.  The  other  ISAR  data  set,  which 
was  used  in  [1],  is  the  6  -  Point  Delta  Wing  experimental  data  set  collected  by  Defence 
Research  and  Development  Canada  -  Ottawa  (DRDC  -  Ottawa)  at  their  radar  range  lo¬ 
cated  at  Shirley’s  Bay,  Ontario.  A  photo  of  the  6  -  Point  Delta  Wing  is  shown  in  Figure 
2  and  an  example  of  its  blurred  and  focused  image  can  be  found  in  Figure  3. 


Unfocused  B727  ISAR  Image 


Cross  Range  -  Doppler  Bins 

Figure  1 .  B727  Simulated  Data  Set  -  Unfocused  and  Focused  Images. 
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Figure  2.  Picture  of  6  -  Point  Delta  Wing  (From  [1].) 
Note  that  one  reflector  is  offset  to  give  a  sense  of  orientation. 


Unfocused  Delta  Wing 


1015  1020  1025  1030  1035  1040 

Cross  Range  -  Doppler  Bins 
Focused  Delta  Wing 


1015  1020  1025  1030 


Cross  Range  -  Doppler  Bins 

Figure  3.  Example  of  an  Unfocused  and  Focused  6  -  Point  Delta  Wing. 

These  two  data  sets  were  used  to  test  the  algorithms  designed  in  this  thesis.  The 
B727  simulated  data  set  is  an  excellent  data  set  since  its  ideal  motion  error  characteristics 
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make  it  perfect  for  testing  the  first  cut  of  an  imaging  algorithm.  The  6  -  Point  Delta 
Wing  data  set  is  an  experimental  data  set  and  it  is  particularly  useful  since  the  entire  data 
set  consists  of  60,000  pulses  in  the  cross-range.  The  60,000  pulses  can  be  cut  into  differ¬ 
ent  size  imaging  intervals  with  each  of  these  intervals  displaying  a  different  amount  of  er¬ 
ror.  Also  the  6  -  Point  Delta  Wing  experimental  data  set  is  the  only  one  of  the  two  data 
sets  that  can  be  used  for  the  3D  motion  detection  portion  of  the  thesis. 

2.  Primary  Research  Questions 

Based  on  the  objectives  above,  there  are  two  research  questions  that  will  be  the 
primary  subject  of  this  thesis: 

a)  Can  the  effective  application  of  a  genetic  algorithm  or  a  particle  swarm 
optimization  algorithm  on  the  parameter  search  of  the  basis  function  for  the  AJTF  algo¬ 
rithm  significantly  reduce  the  computation  burden  required  to  focus  the  ISAR  image? 

b)  Can  the  fast  AJTF  algorithm  found  in  the  first  part  of  this  thesis  be  applied 
to  the  3D  motion  detection  problem  in  order  to  rapidly  sort  the  good  imaging  intervals, 
which  obey  the  2D  motion  assumption,  from  the  imaging  intervals  that  contain  a  high  de¬ 
gree  of  3D  motion  and  cannot  be  focused? 

3.  Organization 

This  thesis  is  organized  into  six  parts. 

i.  Chapter  II  introduces  ISAR  signal  processing  to  the  extent  required  to  un¬ 
derstand  the  Adaptive  Joint  Time-Frequency  Algorithm  (AJTF); 

ii.  Chapter  III  discusses  time-frequency  transforms  and  the  AJTF  algorithm, 

as  found  in  [1],  [2]  and  [3]; 

iii.  Chapter  IV  presents  the  Genetic  Algorithm  (GA)  and  the  Particle  Swarm 
Optimization  (PSO)  Algorithm  that  were  written  to  optimize  the  AJTF  al¬ 
gorithm  along  with  the  ISAR  imaging  results; 

iv.  Chapter  V  presents  the  algorithm  for  3D  motion  detection  and  ISAR  imag¬ 
ing  results  of  imaging  intervals  determined  to  possess  a  high  and  low  de¬ 
gree  of  3D  motion. 
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v.  Chapter  VI  summarizes  the  results  that  were  achieved  in  this  thesis  and  of¬ 
fers  conclusions  and  recommendations  for  future  work;  and 

vi.  The  Appendix  lists  the  MatLab  code  that  was  written  for  this  thesis. 
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II.  INTRODUCTION  TO  INVERSE  SYNTHETIC  APERTURE 
RADAR  (ISAR)  PROCESSING  AND  THE  FOCUSING  PROBLEM 


This  chapter  introduces  the  basic  equations  required  to  understand  ISAR  imaging 
and  the  focusing  problem,  which  occurs  when  non- stationary  signals  such  as  time- 
varying  Doppler  are  present  in  the  ISAR  return  signal.  It  shows  how  Doppler  resolution 
is  related  to  the  coherent  processing  interval  (CPI)  and  how  this  relates  to  the  amount  of 
motion  error  in  an  image. 

For  the  purposes  of  this  thesis,  the  ISAR  return  signal  is  at  the  point  in  the  proc¬ 
essing  where  all  of  the  point  scatterers  have  been  placed  in  their  proper  range  bins  and  a 
fast  Fourier  transform  (FFT)  will  focus  the  scatterer  in  their  appropriate  cross-range  bins. 
The  ISAR  signal  without  error  can  be  represented  as  [1,3]: 

s(x,t)  =  YjAk  exp  j-/— ^-[R0  +  xk  cos +  yk  sin(0(f))]l,  (1) 

k=\  l  C  J 

where  x  is  the  range  bin,  Nk  is  the  number  of  scatterers  in  the  range  bin  x,  Ak  is  the  mag¬ 
nitude  of  the  radar  return  signal,  R0  is  the  distance  to  the  rotational  center  of  the  target 
which  is  constant  during  the  coherent  processing  interval,  (xk,yk)  is  the  down-range  and 

cross-range  distance  from  the  target’s  rotational  center  to  the  kth  point  scatterer.  Since  no 
motion  error  is  being  assumed  at  this  point,  6(t)  can  be  written  as  [1],  [3]: 

e{t)  =  e0+nt,  (2) 

where  Q  is  the  constant  angular  velocity  of  the  radar  target  from  which  the  Doppler 
cross-range  location  of  the  scatterer  is  extracted.  If  6(t)  is  assumed  to  be  small,  Equa¬ 
tion  (1)  reduces  to  [3]: 

s(x, t)  =  YJ  Ak  exp +Xk+  (Q,) yk ]1  (3) 

k=\  l  C  J 

which  can  easily  be  focused  by  the  fast  Fourier  transform  since  the  extracted  Doppler  is 
time-invariant.  The  time-invariance  of  the  Doppler  is  critical  for  successful  focusing  us¬ 
ing  the  fast  Fourier  transform. 
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The  second  equation  of  importance  deals  with  the  Doppler  resolution,  A fD,  and 
its  associated  effect  on  the  cross-range  resolution,  A rcr .  Doppler  resolution  is  associated 
with  the  coherent  processing  interval,  TCPI ,  such  that  [3]: 


4/d  = 


(4) 


which  states  that  as  the  coherent  process  interval  is  increased,  the  Doppler  resolution  in¬ 
creases.  Similarly,  resolution  in  cross-range  can  be  defined  as  [3]: 


A  r 


< 


2Q/0 


¥o  = 


2Q/0 


1  CPI 


(5) 


Therefore,  from  Equation  (5),  it  is  clear  that,  as  the  coherent  processing  interval  is  in¬ 
creased,  the  size  of  the  Doppler  bins  in  the  cross-range  becomes  smaller  and  it  is  possible 
to  resolve  two  different  point  scatterers  in  the  same  range  bin  even  if  they  are  only 
slightly  separated  in  Doppler.  In  this  thesis,  the  number  of  pulses  in  the  cross-range  is 
analogous  to  TCPI .  A  greater  TCPI  leads  to  a  greater  number  of  pulses  in  the  cross-range 

and,  subsequently,  an  ISAR  signal  with  greater  Doppler  separation  between  point  scatter¬ 
ers  that  are  collocated  in  the  same  range  bin. 

If  Equation  (3)  held  for  real  ISAR  returns  from  fast  moving  targets,  it  would  make 
sense  to  use  the  maximum  number  of  radar  pulses  (by  maximizing  the  TCPI )  to  form  the 

high  resolution  image.  However  if  the  entire  60,000  pulses  ( TCPI  =  30  seconds)  of  the  6  - 

Point  Delta  Wing  experimental  data  set  is  focused  using  the  fast  Fourier  transform,  the 
resulting  image,  as  shown  in  Figure  4,  is  extremely  blurred  and  bears  no  resemblance  to 
the  actual  target  in  Figure  2.  Even  if  a  more  reasonable  TCPI  is  chosen,  such  as  a  TCPI  =  1 

second  with  2048  pulses  in  the  cross-range,  the  motion  error  is  still  significant  enough  to 
cause  the  image  to  blur  beyond  recognition.  This  can  be  seen  by  referring  to  the  unfo¬ 
cused  image  in  Figure  3. 

If,  instead,  the  TCPI  is  reduced  to  0.128  seconds  so  that  there  are  now  only  256 

pulses  in  the  cross-range,  motion  error  can  be  avoided.  The  side  effect  of  this  choice  is 
that  the  cross-range  resolution  is  greatly  reduced.  Most  imaging  intervals  will  generate 
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images  that  resemble  those  found  in  Figure  5.  These  images  possess  no  noticeable  error 
in  the  cross-range  but  also  possess  no  resolution  in  the  cross-range  with  most  point  scat¬ 
tered  collapsed  together  on  top  of  each  other  in  the  center  Doppler  bin. 


6  -  Paint  Deta  Wing  Data  Set 
ftfl.mi'i  Pule;**:  in  Or-rw^Eangn-  T  »  gfl  ^rfirwte 


Figure  4.  60,000  Pulse  6  -  Point  Delta  Wing  Data  Set  Focused  with  the  Fast  Fou¬ 

rier  Transform. 


6  -  Point  D&fta  Wing  Data  Set 


m  130 

Cross  Rang*  ■■  Rang*  Bins 


Figure  5.  256  Pulse  Imaging  Interval  6  -  Point  Delta  Wing  Data  Set  Focused  with 

the  Fast  Fourier  Transform. 
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Therefore,  our  motivation  for  a  fast  and  effective  motion  compensation  algorithm 
becomes  obvious.  Any  imaging  interval  with  a  TCPI  that  leads  to  good  separation  in  Dop¬ 
pler  of  the  point  scatterers  within  the  same  range  bin  will  possess  excessive  motion  error, 
time- varying  Doppler  and  a  blurred  image.  An  imaging  interval  with  a  small  TCPI  will 

have  no  error,  time-invariant  Doppler  but  also  no  separation  of  point  scatterers  in  the 
cross-range.  What  is  required  is  an  algorithm  that  can  remove  the  time -varying  Doppler 
from  an  imaging  interval  of  sufficiently  long  TCPI  such  that  the  resulting  image  is  both 
well  focused  and  well  separated  in  Doppler. 

With  our  requirement  for  an  algorithm  to  correct  motion  error  in  an  ISAR  signal 
now  clear,  we  can  proceed  to  Chapter  III  where  time-frequency  transforms  and  the  Adap¬ 
tive  Joint  Time-Frequency  algorithm  are  introduced.  Time-frequency  transforms  and  the 
AJTF  are  both  crucial  tools  in  achieving  effective  motion  compensation  in  an  ISAR  sig¬ 
nal. 
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III.  INTRODUCTION  TO  TIME-FREQUENCY  TRANSFORMS 
AND  THE  ADAPTIVE  JOINT  TIME-FREQUENCY  ALGORITHM 


As  stated  in  Chapter  II,  in  order  to  form  good  ISAR  images  that  are  well  resolved 
in  the  cross-range,  there  is  a  requirement  for  a  coherent  processing  interval,  TCPI,  that  is 

sufficiently  long  such  that  the  cross-range  resolution,  A rcr ,  is  acceptably  small.  Unfortu¬ 
nately  a  TCPI  that  leads  to  a  small  A rcr  also  leads  to  blurring  in  the  cross-range  since  the 

Doppler  of  the  target’s  scatterers  is  no  longer  stationary,  which  is  a  critical  requirement 
of  the  FFT.  In  this  chapter,  time-frequency  transforms  are  introduced  so  that  the  non¬ 
stationary  nature  of  the  time-variant  Doppler  can  be  understood.  Also,  the  AJTF  algo¬ 
rithm,  which  is  an  effective  way  to  compensate  for  time-variant  Doppler  is  introduced. 

A.  TIME-FREQUENCY  TRANSFORMS 

As  stated  in  [1],  [3]  and  [5],  the  one  weakness  of  the  fast  Fourier  transform  (FFT) 
is  that  it  loses  all  information  on  how  the  signal  changes  with  time.  As  explained  in 
Chapter  II,  ISAR  signal  processing  that  uses  the  assumption  of  time-invariant  Doppler 
during  the  coherent  processing  interval,  TCPI ,  will  inevitably  introduce  severe  blurring 

into  the  final  high  resolution  image.  Therefore,  it  becomes  necessary  to  introduce  the 
time-frequency  transform  and  how  it  can  be  used  to  create  an  analogy  between  what  is 
occurring  in  the  time-frequency  plane,  the  Fourier  Spectrum  and  the  radar  image. 

1.  Short  Time  Fourier  Transform 

As  discussed  in  [5],  the  short  time  Fourier  transform  (STFT)  is  a  unique  and 
straightforward  way  of  introducing  a  time  dependency  into  the  Fourier  transform.  This 
is  done  by  pre-windowing  the  FFT  around  a  certain  time  instance  and  can  be  defined  as 
[5]: 

+00 

S,(t,v,h)=  J  s,(u)h\u-t)e-/2’™du  (6) 

—00 

where  Sx  is  the  STFT  of  range  bin  x,  sx  is  the  ISAR  signal  for  range  bin  x  and  h*  (u  - 1)  is 

the  complex  conjugate  of  the  windowing  function,  which  for  purposes  of  this  thesis  is  the 
default  Hamming  window.  Figure  6  shows  the  time-domain  waveform  of  an  arbitrary 
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ISAR  radar  signal  and  the  Hamming  window.  To  generate  the  STFT,  the  Hamming  win¬ 
dow  is  moved  across  the  radar  signal  and  at  each  time  instance  the  FFT  is  taken  of  the  re¬ 
sult  of  the  multiplication  of  the  window  with  the  ISAR  radar  signal. 

Using  the  time-frequency  toolbox  from  [5],  the  STFT  of  an  ISAR  radar  signal  can 
be  generated  and  the  analogy  between  what  is  shown  by  the  STFT,  seen  in  the  power 
spectrum  of  an  ISAR  range  bin  and  displayed  in  the  radar  image,  can  be  illustrated.  Fig¬ 
ure  7  shows  these  three  graphs  when  taken  from  a  blurred  ISAR  image.  By  examining 
range  bin  3 1  in  the  close-up  of  the  B727  blurred  ISAR  image  (left),  it  can  be  seen  that 
there  appears  to  be  two  scatterers  smeared  across  several  Doppler  bins  in  the  cross-range. 
The  same  deduction  can  be  made  by  examining  the  power  spectrum  (upper  right),  as  it 
can  be  seen  that  there  are  two  Doppler  smeared  point  scatterers  located  in  this  range  bin. 
Finally,  the  STFT  (lower  right)  can  be  examined.  From  the  STFT,  it  can  be  very  clearly 
seen  that  there  are  two  point  scatterers  located  in  range  bin  3 1 .  The  STFT  also  shows  the 
nature  of  the  time- varying  Doppler  of  the  two  point  scatterers  as  the  ISAR  signal  is  inte¬ 
grated  from  pulse  to  pulse  during  the  coherent  processing  interval.  The  first  point  scat- 
terer’s  Doppler  decreases  from  pulse  to  pulse  while  the  second  point  scatterer’s  Doppler 
increases  from  pulse  to  pulse. 

Figure  8  shows  the  same  information  as  Figure  7;  however,  it  is  for  a  focused 
ISAR  image.  Examining  the  close-up  ISAR  image  from  a  focused  B727  data  set,  it  can 
be  seen  that  the  point  scatterers  in  range  bin  3 1  have  collapsed  in  the  cross-range  and  are 
well  focused.  It  can  also  be  seen  that  these  scatterers  have  been  shifted  in  Doppler. 

Since  this  shift  is  applied  equally  across  the  image,  it  has  no  adverse  affect.  More  impor¬ 
tant  are  the  results  that  are  seen  in  both  the  STFT  and  the  power  spectrum  of  range  bin 
3 1 .  The  STFT  shows  that  the  two  point  scatterers,  which  are  now  lines  with  zero  slope  in 
the  time-frequency  plane,  no  longer  possessing  time-varying  Doppler  during  the  coherent 
processing  interval.  This  characteristic  is  echoed  in  the  power  spectra  of  the  two  point 
scatterers  which  now  come  to  sharp  points  rather  than  exhibiting  smearing  of  their  power 
spectra  across  several  Doppler  bins  as  occurred  in  the  unfocused  image. 
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Time  Signal  of  Range  Bin  31  and  the  Hamming  Window 


Figure  6.  ISAR  signal,  s^,  of  Range  Bin  31  and  the  Hamming  Window  (After 

[5]-) 
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Figure  7.  Blurred  ISAR  Image  (close-up  of  Range  Bin  31  in  B727  data  set),  Power 

Spectrum  of  Range  Bin  31  and  its  STFT. 
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Focused  Image  -  B727 
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Figure  8.  Focused  ISAR  Image  (close-up  of  Range  Bin  3 1  in  B727  data  set), 
Power  Spectrum  of  Range  Bin  3 1  and  its  STFT. 


In  concluding  this  section,  it  is  important  to  reiterate  the  following  points  as  they 
will  become  critical  in  chapter  4  when  image  focusing  is  discussed: 

a)  The  range  bins  of  a  blurred  ISAR  image  will  have  point  scatterers  that  ap¬ 
pear  as  lines  with  slopes  (or  even  curved  for  higher  order  motion)  when  transformed  to 
the  time-frequency  plane  and  the  slope  is  an  indication  of  the  presence  of  time-varying 
Doppler.  This  same  range  bin  will  have  the  peaks  of  its  power  spectrum,  which  are  pro¬ 
duced  by  each  point  scatterer  in  the  range  bin,  smeared  across  several  Doppler  bins  when 
time-varying  Doppler  is  present  in  the  signal. 

b)  The  range  bins  of  a  focused  ISAR  image  will  have  point  scatterers  that 
appear  as  lines  with  a  slope  of  zero  in  the  time-frequency  plane.  The  power  spectrum  of 
these  same  range  bins  will  have  sharp,  well  defined  peaks  located  solely  in  their  appro¬ 
priate  Doppler  bins  with  minimum  spillage  into  adjacent  bins. 


14 


2.  Weaknesses  of  the  STFT  and  Time  -  Frequency  Transform 

Despite  the  usefulness  of  the  STFT  and  other  time  -  frequency  transforms,  as 
listed  in  [5]  and  [3],  there  are  two  primary  weaknesses  which  will  be  of  concern  in  the  fo¬ 
cusing  problem  of  this  thesis.  The  first  weakness  is  the  computational  time  of  the  STFT, 
which  is  the  fastest  time-frequency  transform  available  in  the  time-frequency  toolbox 
from  [5],  Calculation  time  for  the  STFT  is  significantly  longer  than  the  FFT  and  this  be¬ 
comes  critical  if  the  STFT  must  be  computed  multiple  times.  The  second  weakness  of 
the  STFT,  as  stated  in  [2]  and  [5],  is  that  it  possesses  poor  resolution  and  cannot  be  si¬ 
multaneously  well  resolved  in  both  time  and  frequency.  This  is  not  a  significant  concern 
in  the  B727  simulated  data  set  as  all  point  scatterers  in  each  range  bin  are  well  separated 
in  Doppler.  However,  for  the  6  -  Point  Delta  Wing  experimental  data  set,  the  point  scat¬ 
terers  are  not  well  separated  in  Doppler  and  it  is  difficult  to  see  how  the  Doppler  of  each 
point  scatterer  changes  with  time.  To  increase  time-frequency  resolution,  as  recom¬ 
mended  in  [3]  and  [5],  a  different  time-frequency  transform  could  be  chosen;  however, 
this  typically  results  in  an  increase  in  computational  time,  or  in  the  case  of  the  Wigner- 
Ville  Distribution  (WVD),  cross-terms  appear  in  the  time-frequency  image.  As  shown  in 
Figure  9,  which  is  a  WVD  of  range  bin  31  of  the  unfocused  B727  simulated  data  set, 
these  cross-terms  can  make  the  interpretation  of  the  time-frequency  image  very  difficult 
and  the  automation  of  line  detection  close  to  impossible. 
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Figure  9.  Wigner-Ville  Distribution  of  Unfocused  Range  Bin  31  (After  [3]  and 

[5]-) 
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B.  ADAPTIVE  JOINT  TIME-FREQUENCY  ALGORITHM 

The  AJTF  algorithm  can  be  found  in  [1],  [2]  and  [3]  and  is  capable  of  removing 
time-varying  Doppler  in  a  radar  signal  that  is  due  to  translational  and  rotational  motion 
error.  As  explained  in  the  references,  the  capability  of  the  AJTF  is  limited  by  the  hold¬ 
ing  of  three  primary  assumptions.  The  first  assumption  is  that  the  ISAR  signal  has  been 
error  corrected  in  the  down-range  direction  and  that  each  point  scatterer  is  in  its  proper 
range  bin.  The  second  assumption  is  that  of  a  rigid  body,  which  implies  that  motion  error 
that  exists  in  one  point  scatterer  will  exist  in  all  point  scatterers  of  the  target.  This  as¬ 
sumption  is  violated  in  the  presence  of  jet  engine  modulation  (JEM),  helicopter  rotors  and 
the  micro-Doppler  phenomenon.  The  third  assumption  is  that  motion  occurs  in  a  2D 
plane  with  respect  to  the  radar.  The  AJTF  algorithm  is  currently  unable  to  correct  for  3D 
motion.  For  this  thesis,  the  first  two  assumptions  will  always  hold.  The  assumption  of 
2D  motion  does  not  hold  for  all  imaging  intervals  of  the  6  -  Point  Delta  Wing  experimen¬ 
tal  data  set  and,  as  a  result  of  the  3D  motion  error,  these  imaging  intervals  will  not  image 
well.  As  stated  in  the  introduction,  it  is  the  second  part  of  this  thesis  that  deals  with  3D 
motion  error  so,  for  now,  it  will  be  assumed  that  all  three  assumptions  hold. 

1.  The  Motion  Model  for  a  Moving  ISAR  Target 

The  motion  model  for  an  ISAR  target  is  introduced  in  [2]  and  further  explained  in 
both  [3]  and  [1],  The  following  explanation  is  drawn  from  all  three  sources  but  primarily 
[1]  and  [3],  Figure  10  shows  the  motion  geometry  of  an  ISAR  target.  Coordinates 
(U,V,W)  is  the  radar  coordinate  system,  while  coordinates  (x,y,z)  is  the  coordinate 

system  used  to  define  the  target  and  (X,Y,Z)  is  the  coordinate  system  translated  from 
the  radar  and  is  used  to  describe  the  target’s  rotation.  It  is  centered  at  the  geometric  cen¬ 
ter  of  the  target,  about  which  the  target  will  rotate.  If  P(x,y )  is  the  location  of  a  point 

scatterer,  then  the  initial  distance  from  the  radar  to  the  klh  point  scatterer  located  on  the 
target  at  t  =  0  is  [3] : 

Rk  =Ro  +  xpcos(0o-a)  +  ypsm(0o-a),  (7) 
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where  R0  is  the  initial  line-of-sight  distance  from  the  radar  to  the  target’s  geometric  cen¬ 
ter,  0O  is  the  initial  angle  between  (X,Y,Z)  and  the  (x,  y,z)  coordinate  system  and  a  is 
the  azimuth  angle  of  the  scatterer  off  of  the  (U,  V,W)  axis. 


Q  rotation 


Figure  10.  Geometry  of  the  Radar  Imaging  of  a  Target  (From  [3].) 

Now,  if  the  target  possesses  a  rotational  motion,  Q ,  about  the  Z  axis  and  a  trans¬ 
lational  motion  with  radial  velocity,  VR ,  then  both  the  range  and  rotational  angle  become 

a  function  of  time  and  Equation  (7)  become  time  dependent.  The  time  dependent 
i?(^)can  be  expanded  into  a  Taylor  series  polynomial  as  [1]: 

R(t)=Ra+v„t+±v't2+±iv;t1+...,  (8) 

dV 

where  VL  =  — —  =  aK  which  is  the  acceleration  of  the  scatterer  in  the  radial  direction  and 

*  dt  R 

the  additional  derivatives  of  VR  are  terms  of  higher  order  translational  motion.  The  time 
dependent  0{t )  can  be  expanded  into  a  Taylor  series  polynomial  as  [2]: 

0(O  =  6>o+Qt  +  ^Q't2+^Q"t3  +  ..„  (9) 
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where  Q'  =  which  is  the  initial  angular  acceleration  of  the  point  scatter  and  the  addi¬ 
tional  derivatives  of  Q  are  terms  of  higher  order  rotational  motion.  Taking  Equations  (8) 
and  (9),  substituting  them  into  Equation  (7)  and  arbitrarily  setting  a  and  0O  to  zero,  we 
have  [2]: 


^*(0  --^o  +  +  {2  +  —  +  xkcos 


f 


1  '-\i. 2  .  1 


A 


Qt  +  — Q'r  +  — Q"f  + ... 
v  2!  3!  j 


+yk  sin 


1  0^2  .  1  rv'*3 


Qt  +  — Q'r  +— Q"f  +... 

V  2!  3! 


(10) 


Now  taking  Equation  (10)  and  calculating  the  Taylor  series  expansion  of  the  sine  and  the 
cosine  about  t  yields  [1]: 


Rk (t)-{R0+xk)+ (Qf*  +  VR)t+^{K  - Q2 xk + ^Nc)/2  +  ••• 


(11) 


which  is  the  time-varying  range  to  a  one  point  scatterer  on  the  target. 

2.  The  ISAR  Signal  and  the  Motion  Compensation  Algorithm 
Now  Equation  (11)  can  be  substituted  into  Equation  (1)  which  gives  the  equation 
for  an  ISAR  return  signal  from  a  collection  of  Nk  point  scatterers  located  in  range  bin  x 
where  the  ISAR  return  signal  is  not  an  ideal,  errorless  signal  [1]: 


.v(x,?)  =  X4-exP]- 


:4^/o 


k-\ 


(R0  +  xk)  +  {O.yk+VR)t  +—(VR-Cl2xk+Q.yk)t2+... 


(12) 


If,  as  in  [1],  the  following  substitutions  are  made: 

a  o  =  +  Xk 

«i=Q yk+vR  (i3) 

a2——(VR  —  Q  xk  +  Qjj.), 


the  following  observations  are  made.  First,  a0  defines  the  coordinate  position  of  the  scat¬ 
terer  on  the  target.  The  second  term,  ax ,  is  the  radar  line-of-sight  translational  motion 

and  constant  rotational  motion  and  is  linearly  dependant  on  time.  If  all  other  terms  are 
negligible,  then  a  simple  fast  Fourier  transform  is  sufficient  to  form  a  focused  image. 
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The  second  term,  a2  (as  well  as  subsequent  terms  if  higher  order  motion  is  present),  con¬ 
tains  V'R  which  is  the  translational  motion  error  and  Q 'yk  which  is  the  rotational  motion 

error  that  has  a  dependency  on  the  cross-range.  The  (Q t)2  can  be  neglected  since  [1] 
states  that  for  centimeter  and  millimeter  wave  radars,  this  term  is  generally  much  less 
than  1. 

Therefore,  we  can  define  the  phases  of  the  uncompensated  range  bin  as  the  phase 
function  [1]: 

O  (^)  =  Uq  +  u^t  +  u2t  +...  (14) 


and  since  the  phases  in  the  range  bin  are  time-varying,  the  instantaneous  frequency  can 
be  expressed  as  [1]: 


m= 


1  ao(r) 
2  n  dt 


(15) 


which  is  also  time-varying.  Therefore,  the  Doppler  frequencies  of  the  point  scatterers  in 
the  range  bin  are  time-varying  and  the  image  is  blurred  in  the  cross-range. 


The  next  step  is  to  focus  on  a  range  bin  and  to  find  the  translational  motion  com¬ 
pensation  basis  function  of  the  form  [3]: 


hPT(t)  =  exp{-j2^-®T(t)} 
=  exp + 


(16) 


where  &T(t )  is  the  phase  function  for  translational  motion  compensation  and  the  parame¬ 
ters  { /n,/i2,/13,...}  are  the  AJTF  search  parameters  for  translational  motion  compensa¬ 
tion.  A  suitable  basis  function  will  closely  resemble  the  prominent  point  scatterer  in  the 
time-frequency  plane.  Figure  1 1  shows  an  STFT  of  the  unfocused  range  bin  from  the 
B727  simulated  data  set  and  the  STFT  of  the  best  basis  function  for  translational  motion 
compensation  (described  later  in  this  chapter).  Once  this  basis  function  has  been  found, 
translational  motion  error  can  be  corrected  by  applying  the  following  equation  from  [1], 
[2]  and  [3]: 
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s(x,?)r  =s(*,f).*conj{Aft-  0)} 


(17) 


where  s(x,t)T  is  the  ISAR  signal  for  all  range  bins  with  translational  motion  error  re¬ 
moved,  s(x,t)  is  the  ISAR  signal  for  all  radar  bins,  .*  defines  array  multiplication  and 

conj  is  the  MatLab  complex  conjugate  operator.  Figure  12  shows  the  same  range  bin  in 
the  time-frequency  plane  after  translational  motion  compensation.  Note  that  the  top  line, 
representing  the  prominent  scatterer  in  the  range  bin  has  had  its  slope  zeroed  which  indi¬ 
cates  the  removal  of  time-varying  Doppler.  The  remaining  signal  is  now  [1]: 

s(x,t)T  =  ZAexp|-7^^[Ax/t+Ay(t6>(0]|  (18) 

where  0{t)  is  as  defined  in  Equation  (9)  above. 
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Figure  1 1 .  Unfocused  Range  Bin  (top)  and  the  Best  Basis  Function  hpT  ( t )  for 

Translational  Motion  Compensation. 
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Range  Bin  31  after  Translational  Motion  Compensation 
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Figure  12.  Time-Frequency  Transform  of  Range  Bin  After  Translational  Motion 

Compensation. 

Now  it  is  required  that  rotational  motion  compensation  be  completed  on  the  sig¬ 
nal.  This  is  done  by  choosing  another  range  bin  with  a  prominent  point  scatterer.  If  there 
is  more  than  one  scatterer  in  a  range  bin,  it  is  possible  to  use  the  same  range  bin;  how¬ 
ever,  the  prominent  point  scatterer  used  for  translational  motion  must  first  be  removed 
from  the  signal.  The  algorithms  written  for  this  thesis  used  separate  range  bins  for  trans¬ 
lational  and  rotational  motion  compensation.  The  first  step  is  to  again  find  a  suitable  ba¬ 
sis  function  as  in  Equation  (16)  and  as  shown  in  Figure  9: 

hPR(t)  =  exp{-j2x-®R(t)}.  (19) 

However,  for  rotational  motion  compensation,  it  is  the  phases  of  the  basis  function  in 
Equation  (19),  &R  (t) ,  that  are  required  [3]: 

=  fi\t +  —fiit2  +  + (20) 
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where  &r(t)  is  the  phase  function  for  rotational  motion  compensation  and  the  parameters 

{/2i,/22,/23,...}  are  the  search  parameters  for  rotational  motion  compensation.  The  phase 

function  of  Equation  (20)  and  the  rotational  angle  in  Equation  (18)  gives  the  relationship 
between  rotational  angle  and  dwell  time  such  that  [3]: 

0(O«©*(O.  (21) 

The  linear  interpolation  function  in  MatLab  can  reformat  the  ISAR  data  so  that  the  radar 
signal  is  uniformly  sampled  in  6 ,  only  linearly  dependent  on  time  and  the  quadratic 
phase  error  is  removed.  The  final  focused  ISAR  signal,  as  given  in  [1],  becomes: 

s(.x>0R  =  iAk  exp|-y^^[Ar/(  +  Ayt0(O]  j  (22) 

where  s(x,t')R  is  the  ISAR  signal  after  translational  and  rotational  motion  compensation, 
6(t')  is  the  uniformly  sampled  rotational  motion  term  that  only  depends  linearly  on  t'  and 
t'  is  the  new  time  variable. 

3.  Selection  of  an  Appropriate  Basis  Function  or  Phase  Function 

In  the  last  section,  it  was  stated  that  suitable  basis  functions,  h„  (t)  and  h„  (t), 

must  be  found  in  order  to  perform  translational  and  rotational  motion  compensation. 

There  are  three  methods  that  can  be  used  to  determine  if  a  basis  function  is  suitable  for 
compensating  for  motion  error.  The  most  common  method  is  the  projection  method 
which  it  is  used  in  [1],  [2],  [3]  and  [6]  as  the  way  to  determine  the  suitability  of  the  basis 
functions.  The  second  method  follows  the  work  in  [16]  and  [17]  by  using  the  STFT  and 
the  Radon  transform  to  determine  basis  function  suitability  by  automatically  recognizing 
the  conditions  of  focus.  The  third  method  of  determining  basis  function  suitability  uses 
the  FFT  to  automatically  recognize  focus  through  the  characteristics  seen  in  the  power 
spectrum.  Note  that  for  these  three  methods,  first  the  best  translational  motion  compen¬ 
sation  parameters,  {/n,/12,/13,...},  are  found  and  translational  motion  compensation  is 
performed  on  the  signal.  After  this  has  been  completed,  the  rotational  motion  compensa¬ 
tion  parameters,  {/21,/22,/23,...},  can  be  found. 
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a.  Projection  Method  to  Determine  hpi  (?)  and  hpR  (?) 

This  approach  is  described  mathematically  in  [1]  as: 

if\\J\ 2,/i3,-}  =  argmax  (s(x,t)hpK(t))  (23) 

where  (s(x,  t)hpx  (t)j  is  the  projection  of  the  basis  function  onto  the  signal  of  range  bin  x 
in  the  time  domain  and  can  be  defined,  as  in  [1],  as: 

(*(*, o*„«)|= «<*•  (24) 

Rotation  motion  compensation  parameters  are  found  similarly  as: 

{fivfnJw)  =  argmax|(s(x,?)r/ipr(?))|  (25) 

where  (s{x,t)ThpTit)^  is  defined  in  a  similar  fashion  as  Equation  (24).  Therefore,  from 

the  above  equations,  the  most  suitable  basis  function  will  be  the  one  with  the  parameters 
{/ii,/i2,/i3,...}  or  {/2i,/22,/23,...}  that,  after  vector  multiplication  between  the  signal  and 

the  basis  function  and  after  taking  the  magnitude  of  the  result,  will  produce  the  greatest 
value.  Though  this  method  is  a  very  popular  method  for  determining  basis  function  suit¬ 
ability,  it  was  observed  that  the  projection  method  did  not  always  lead  to  an  adequate  ba¬ 
sis  function.  To  address  this,  the  next  two  methods  are  discussed. 

b.  Automatic  Recognition  of  Focus  using  the  STFT  and  the  Radon 
Transform 

The  Radon  transform  has  been  previously  used  in  1SAR  image  formation 
in  [16]  and  [17].  In  this  thesis,  the  Radon  transform  was  used  to  determine  the  suitability 
of  the  basis  and  phase  functions,  hfT  (/)  and  0/;(?) ,  and  is  a  replacement  for  the  previ¬ 
ously  explained  projection  method.  Information  on  the  calculation  of  the  Radon  trans¬ 
form  is  found  in  [13]  and  detection  of  lines  in  the  time-frequency  plane  is  introduced  in 
[5],  This  method  finds  its  validity  in  the  argument,  as  presented  in  Chapter  III,  that  a  fo¬ 
cused  scatterer  which  has  had  its  motion  error  and,  consequently,  its  time -varying  Dop¬ 
pler  removed,  will  be  a  line  of  zero  slope  in  the  time-frequency  plane.  The  first  step  of 
this  method  is  to  take  the  current  search  parameters  { /11,/12,/13,...}  and  formulate  the  ba- 
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sis  function,  h  (t)-  (If  translational  motion  compensation  has  already  been  completed, 
the  search  parameters  {/2i,/22,/23,...}  are  used  to  form  the  phase  function,  Q/t  (/),  for  ro¬ 
tational  motion  compensation.)  Once  hfR  (t)  or  (?)  has  been  formulated,  motion  com¬ 
pensation,  either  translational  or  rotational  as  applicable,  is  performed  on  the  range  bin. 

It  is  only  necessary  to  perform  motion  compensation  on  the  range  bin  under  examination 
during  the  basis  and  phase  function  searches  since  the  rigid  body  assumption  states  that 
motion  error  is  constant  from  range  bin  to  range  bin.  The  next  step  is  to  take  the  STFT  of 
the  resulting  range  bin  and  to  calculate  the  Radon  transform  of  the  resulting  image. 

The  Radon  transform  is  a  proven  method  for  detecting  lines  in  images  (see 
[5]  and  [13]).  It  is  able  to  return  parameterized  information  on  the  lines  found  in  the  im¬ 
age.  Figure  13  is  the  plot  of  the  Radon  transform  of  the  STFT  of  range  bin  3 1  in  the  B727 
simulated  data  set  in  its  unfocused  state,  after  translational  motion  compensation  and  af¬ 
ter  rotational  motion  compensation.  The  images  of  the  time-frequency  transforms  have 
already  been  displayed  for  these  three  states  in  Figures  11,12  and  8,  respectively.  For 
the  purposes  of  this  method,  only  the  6  parameter  of  the  Radon  transform  is  important  as 
a  value  of  6  =  90°  indicates  a  flat  line  in  the  subject  time-frequency  plane.  Figure  13 
demonstrates  that,  as  lines  in  the  time-frequency  plane  are  flattened  to  zero  slope  when 
correct  values  of  {/11,/12,y]3,...}  and  {/2i,/22,/23,...}  are  found  and  motion  compensation 

is  correctly  applied,  the  peaks  in  the  Radon  transform  graphs  lie  on  the  6  =  90°  line.  The 
final  Radon  transform,  which  is  taken  from  the  STFT  of  a  focused  image  shows  that  there 
are  two  lines  in  the  image,  one  laying  on  the  0  =  90°  line  and  the  other  on  the  0  =  89° 
line.  Therefore,  detecting  focus  is  simply  a  matter  of  searching  for  the  parameters 
{/ii,/i2,/i3,...}  and  {/21,/22,/23,...}  that  cause  peaks  along  or  very  near  the  6  =  90°  line 
when  motion  compensation  is  applied. 

The  one  significant  drawback  of  this  method  which  unfortunately  makes  it 
currently  unusable,  is  that  both  the  STFT  and  the  Radon  transforms  take  too  long  to  cal¬ 
culate  and  even  the  most  efficient  of  search  algorithms  must  calculate  these  transforms 
many  times.  Therefore,  this  method  cannot  be  currently  used  in  ISAR  image  focusing  in 
a  real-time  application.  However,  with  the  ever  increasing  computational  power  in  to- 
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day’s  CPUs,  it  may  one  day  be  a  suitable  for  determining  basis  function  suitability.  For¬ 
tunately,  the  research  into  this  area  led  to  the  realization  of  a  method  that  uses  the  FFT. 


Radon  Transform  o(  Uncompensated  Range  Bin  *  10*  Radon  Transform  oS  Transitional  Motion  Compensated  Range  Bin  x  10* 
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Figure  13.  Radon  Transform  of  Uncompensated  and  Compensated  Range  Bins. 


c.  Automatic  Recognition  of  Focus  using  the  FFT 
This  method  was  developed  soon  after  it  was  realized  that  the  method  us¬ 
ing  the  STFT  and  the  Radon  transform  could  not  be  practically  implemented.  The  FFT 
method  finds  its  validity  in  the  explanation  given  in  Chapter  III.A  in  that  a  time- 
frequency  transform  of  a  range  bin  with  motion  error  will  have  point  scatterers  mapped 
out  as  tilted  (or  possibly  curved)  lines  and  a  power  spectrum  which  has  areas  of  maxi¬ 
mum  value  smeared  across  several  Doppler  bins.  A  focused  range  bin  will  have  a  scat¬ 
terers  mapped  out  as  lines  with  zero  slope  in  the  time-frequency  plane  and  a  power  spec¬ 
trum  where  maximums  are  sharp  points  in  there  respective  Doppler  (i.e.,  cross-range) 
bins  (refer  to  Figures  7  and  8  for  images  that  illustrate  this). 

The  FFT  method  starts  the  same  as  the  previous  method  in  Subsection 
3.B.b  of  this  Chapter,  in  that  the  range  bin  under  examination  is  subjected  either  to  trans- 
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lational  or  rotational  compensation  using  the  test  parameters  {/n,/12,/13,...}  or 
{/2i,/22,/23,...},  as  applicable.  The  power  spectrum  of  the  range  bin  x  is  then  calculated 
using  the  familiar  equations: 


and 


FFTfs(i)rL)|2 


(26) 


(27) 


where  5  ( t)T  \  is  the  ISAR  signal  of  range  bin  x  after  translational  motion  compensation 
and  S  (/)r|  v  is  its  power  spectrum,  .v  (V ) w  |  r  is  the  ISAR  signal  of  range  bin  y  after  rota¬ 
tional  motion  compensation  and  *S'(/)/j|  is  its  power  spectrum.  Now  to  determine  the 

suitability  of  the  basis  function,  hpr  (j) ,  or  the  phase  function,  &R  (t) ,  it  is  necessary  to 
determine  the  following: 


and 
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Equations  (28)  and  (29)  state  that  the  best  search  parameters  would  be  ones  that  place  the 
maximum  percentage  of  the  power  spectrum  found  in  the  entire  range  bin  in  a  single 
Doppler  bin.  A  simple  maximum  search  will  not  work  since  incorrect  values  of  the 
search  parameters  can  generate  a  high  peak  in  a  Doppler  bin  which  is  much  greater  than 
what  will  be  found  if  the  image  were  focused,  but  there  will  also  be  a  very  significant 
amount  of  the  power  spectrum  in  adjacent  bins  and  accordingly  the  image  will  be  ex¬ 
tremely  blurred.  Taking  the  average  is  very  effective  at  avoiding  this  situation  because 
the  averaging  process  will  score  these  sets  of  parameters  very  low.  An  example  of  a 
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range  bin,  motion  compensated  with  two  different  sets  of  parameters,  one  that  does  not 
focus  the  image  and  the  other  that  focuses  the  image,  is  shown  in  Figure  14.  This  auto¬ 
matic  recognition  of  focus  technique  using  the  FFT  is  very  fast  and  is  the  technique  that 
is  used  in  the  search  algorithms  in  Chapter  IV  for  determining  the  suitability  of  a  basis 
function’s  or  phase  function’s  search  parameters. 

Power  Spectral  Density 
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Figure  14.  Example  of  a  Unfocused  and  Focused  Range  Bin  and  Their  Score. 

C.  CHAPTER  SUMMARY 

In  this  chapter,  two  key  concepts,  which  are  critical  to  Chapter  IV  and  the  overall 
problem  of  ISAR  focusing,  were  introduced.  First,  the  mathematical  model  of  the  AJTF 
algorithm  was  presented  and  it  is  this  algorithm  that  enables  time-varying  Doppler  due  to 
target  motion  to  be  removed  from  the  ISAR  signal.  Secondly,  the  FFT  method  of  deter¬ 
mining  basis  and  phase  function  suitability  was  presented.  Chapter  IV  uses  this  method 
as  the  core  of  the  GA  and  the  PSO  search  algorithms. 


Power  Spectral  Density 
xIO18  of  Focused  Range  Bin 
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IV.  GENETIC  ALGORITHM  AND  PARTICLE  SWARM 
OPTIMIZATION  AS  SEARCH  ALGORITHMS 


This  chapter  covers  the  Genetic  Algorithm  (GA)  and  the  Particle  Swarm  Optimi¬ 
zation  (PSO)  algorithm  and  their  application  to  extracting  translational  and  rotational  mo¬ 
tion  correction  parameters  from  the  solution  space.  The  chapter  concludes  by  offering  a 
comparison  between  the  two  algorithms. 

A.  INTRODUCTION  TO  THE  SOLUTION  SPACE  AND  THE  EXHAUSTIVE 
SEARCH 

Before  beginning  the  discussion  on  the  GA  and  the  PSO  algorithms  it  is  important 
to  briefly  discuss  the  construct  of  the  solution  space  which  must  be  searched  and  the  ex¬ 
haustive  search  method,  as  described  in  [1]  and  [3],  that  is  typically  used  to  find  the 
search  parameters  or  {fivfn’fn’-]-  In  the  AJTF  algorithm  for  ISAR 

image  focusing,  the  parameters  fn  and  f2l  are  always  found  using  the  FFT  such  that: 

fu  =  DopplerBin  jarg  max  |  FFT  (/I)|JV  )|2)|  (30) 

and 

/21  =  DopplerBin  jarg  max  |  FFT  (31) 

which  means  that  fn  and  f2l  are  always  equated  to  the  Doppler  bin  possessing  the  maxi¬ 
mum  power  spectral  density  in  range  bins  x  and  v.  The  remaining  parameters, 

{/i2,/i3,...}  and  {/22,/23,...},  must  be  searched  and  each  parameter  can  possess  any  real 

valued  number  lying  between  -N  to  +N,  where  N  is  the  number  of  pulses  used  to  form  the 
ISAR  image  in  cross-range.  Therefore,  as  the  number  of  pulses  that  are  used  to  form  an 
image  increases,  so  does  the  solution  space  and  the  computational  intensity  required  to 
achieve  focus. 

Figures  15  and  16  shows  a  map  of  the  solution  space  that  must  be  searched  for  the 
B727  simulated  data  set  when  the  search  is  limited  to  a  two-dimensional  search  space 
(i.e.,  3rd-degree  motion  error  which  requires  a  search  for  parameters  {/12,/22}  and 

{/22,/23}).  For  translational  motion  compensation,  it  can  be  seen  that  the  solution  space 


29 


is  fractured,  with  many  areas  of  separate  local  maximums.  Because  of  this  fracturing, 
conventional  searches,  such  as  a  gradient  search,  may  have  problems  finding  the  global 
maximum  and  instead  converge  to  a  local  maximum  [9,14],  This  local  maximum  may 
not  be  an  acceptable  set  of  parameters  for  focusing  the  image.  The  rotational  motion  so¬ 
lution  space  is  not  nearly  as  fractured  but  it  does  have  a  couple  of  valleys  which  may 
cause  a  gradient  search  to  mistakenly  converge  to  a  local  maximum. 

Solution  Space  for  B727  -  Translational  Motion  Compensation 


Figure  15.  2D  and  3D  Visualization  of  B727  Solution  Space  for  Translational  Mo¬ 
tion  Compensation. 
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Solution  Space  for  B727  -  Rotational  Motion  Compensation 
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Figure  16.  2D  and  3D  Visualization  of  B727  Solution  Space  for  Rotational  Motion 

Compensation. 

A  common  method  for  searching  this  solution  space,  as  detailed  in  [1],  [2]  and 
[3],  is  the  exhaustive  search.  The  exhaustive  search  is  typically  done  by  finding  fn  as  in 
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Equation  (30)  and  then  stepping  through  fn  for  the  span  of  the  solution  space  (while  the 
additional  parameters,  /13,/14  and  higher  terms  are  set  to  zero)  and  finding  the  best  term 

that  satisfies  Equation  (28).  Once  this  term  is  found,  the  next  term  is  stepped  through  the 
solution  space,  with  previous  terms  set  to  their  optimal  values,  until  again  the  maximum 
of  Equation  (28)  is  found.  This  is  repeated,  as  required,  for  all  of  the  search  parameters 
until  the  best  basis  function,  hfT  (t) ,  is  found  and  translational  motion  compensation  is 

performed  on  the  ISAR  signal.  This  identical  methodology  is  then  repeated  to  find  the 
best  phase  function,  0/(  (/). 

Though  this  is  an  effective  method  of  finding  the  parameters  of  motion  compensa¬ 
tion,  it  suffers  from  two  drawbacks.  The  first  is  that  when  searching  for  the  third-order 
parameters  (i.e.,  /13  or  /23  and  above),  the  exhaustive  search  may  not  find  a  true  global 

maximum  since  it  is  doing  a  ID  search  of  a  parameter  with  all  other  parameters  fixed,  de¬ 
spite  the  fact  that  the  parameters  are  not  orthogonal.  Secondly,  even  with  this  reduced 
search  space,  it  is  still  very  computationally  intensive  and  time  consuming  to  step  though 
all  possible  answers.  One  way  of  measuring  computational  intensity  is  to  determine  the 
number  of  times  the  cost  function  must  be  calculated.  In  terms  of  this  exhaustive  search 
algorithm,  the  cost  function  is  Equation  (28)  and  for  the  exhaustive  search  it  must  be  cal¬ 
culated  Ncosl_function  times  as  determined: 

^ cost  _  function  ^ cross  _  range  f  ^ (fn  ,/j3 ^ (f 22  s/21  >•■■)  ’  0^) 

where  Ncross  mnge  is  the  number  of  pulses  in  the  cross-range,  Nl  f^  ^  ,  is  the  number  of 

search  parameters  for  translational  motion  compensation  and  N(  f^  ^  .  is  the  number  of 

search  parameters  for  rotational  motion  compensation.  Therefore,  for  the  B727  simulated 
data  set  with  256  pulses  in  the  cross-range  and  3rd-degree  motion  compensation  (required 
to  find  two  search  parameters  for  both  translational  and  rotational  motion  compensation), 
the  number  of  cost  function  calculations  is: 

^ cost  _  function  =  (2  •  256)  •  2  •  2  =  2048.  (33) 
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Therefore,  when  the  efficiency  of  the  GA  and  PSO  algorithm  is  measured,  it  will 
be  evaluated  not  only  by  its  run  time  but  also  by  the  number  of  cost  function  calculations 
it  requires  to  converge.  Minimizing  the  number  of  cost  function  calculations  is  the  key  to 
an  efficient  algorithm  and  a  clear  measure  of  how  much  better  the  algorithm  is  than  the 
simple  exhaustive  search. 

B.  GENETIC  ALGORITHM  PARAMETER  SEARCH 

Genetic  algorithms,  first  developed  in  1989  by  Goldberg,  whose  work  is  pub¬ 
lished  in  [9],  are  particularly  well  suited  to  searching  a  large  solution  space  that  contains 
one  global  maximum  buried  amid  many  local  minima  and  which  may  be  discontinuous. 

A  GA  search  attempts  to  model  natural  behavior  where  strong  traits  will  propagate 
through  a  breeding  population  while  weak  traits  will  die  off  and  eventually  be  removed 
from  a  population  after  a  number  of  generations.  One  outstanding  feature  of  a  well- 
written  GA  is  that  it  will  rapidly  converge  very  close  to  the  global  maximum.  The  disad¬ 
vantage  is  that  once  it  has  converged  to  this  point,  typically  only  small  improvements  will 
be  made  to  the  found  global  maximum  and  this  only  after  many  generations.  This  trait  is 
illustrated  in  Figure  17  and  is  a  typical  feature  of  any  evolutionary  search.  The  GA  ex¬ 
periences  a  rapid  increase  in  the  best  score  of  its  found  global  maximum  in  the  first  23 
generations  of  the  evolutionary  search.  The  next  32  generations  are  characterized  by  sig¬ 
nificant  but  smaller  increases  to  the  score  of  the  found  global  maximum.  The  remaining 
generations  are  characterized  by  only  a  minimal  and  hardly  noticeable  increase  to  the 
score  of  the  found  global  maximum.  Therefore,  using  this  argument,  it  is  obvious  to 
state  that  with  GAs,  convergence  to  the  absolute  and  unique  maximum  in  not  guaranteed. 
However,  the  probability  that  convergence  has  occurred  increases  with  the  user’s  will¬ 
ingness  to  wait  for  the  algorithm  to  process  successive  generations.  Looking  again  at 
Figure  17,  it  can  be  stated  with  confidence  that  after  the  70th  generation,  the  GA  has  con¬ 
verged  either  at  or  near  the  true  global  maximum. 

In  ISAR  imaging,  there  is  no  knowledge  of  the  location  of  the  true  global  maxi¬ 
mum  nor  is  there  any  way  to  predict  the  likeliness  of  finding  a  better  global  maximum. 

As  successive  generations  occur,  the  probability  of  finding  a  better  global  maximum  di¬ 
minishes.  Also,  since  the  goal  is  to  rapidly  find  acceptable  search  parameters  that  will 
focus  the  ISAR  image,  exiting  the  search  as  expediently  as  possible  becomes  critical. 
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Therefore,  the  rapid  convergence  (i.e.,  increase  in  the  score  of  the  found  global  maxi¬ 
mum)  of  the  GA  early  on  in  the  search  can  be  exploited.  Combining  this  rapid  conver¬ 
gence  with  the  ability  to  optimally  decide  when  to  exit  the  search  will  allow  for  rapid  fo¬ 
cusing  of  ISAR  images.  The  work  of  Li  and  Ling  in  [6]  is  an  excellent  first  publication 
of  a  GA  applied  to  ISAR  focusing  and  was  used  as  a  reference  in  the  writing  of  the  GA 
code  for  this  thesis.  In  particular,  [6]  states  that  real-valued  GAs  significantly  outperform 
binary  coded  GAs.  Because  of  this,  the  GA  for  this  thesis  will  forego  exploring  binary- 
coded  GAs  and  investigate  only  a  real-valued  GA.  Another  reference  source  was  [10], 
which  provides  a  good  explanation  of  the  flow  and  construct  of  a  GA. 


Genetic  Algorithm 

Best  Score  and  Average  Score  vs  Generation 


Figure  17.  Illustration  of  a  GA’s  Rapid  Convergence. 

1.  Sequence  of  the  Genetic  Algorithm 

The  genetic  algorithm  follows  a  set  sequence  that  repeats  until  done  (see  Figure 
18)  with  each  cycle  referred  to  as  a  generation.  The  definition  of  each  step  is  as  follows: 
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Figure  18.  Flow  Chart  of  the  ISAR  Genetic  Algorithm  (After  [6].) 

a.  Define  the  Population 

The  first  step  of  the  GA  is  to  generate  a  population  of  suitable  size  to  en¬ 
sure  sufficient  randomness  in  the  initial  group  of  parents.  It  is  important,  however,  to  en¬ 
sure  that  the  population  is  not  so  large  that  it  becomes  unwieldy  and  slows  down  the  algo¬ 
rithm  or  that,  after  only  a  few  generations  of  the  GA,  the  cost  function  has  been  calcu¬ 
lated  more  times  than  it  would  have  been  if  an  exhaustive  search  had  been  used.  How¬ 
ever,  in  the  ISAR  case,  the  ideal  population  size  will  be  function  of  the  number  of  pulses 
used  in  the  cross-range,  the  number  of  search  parameters  due  to  the  order  of  motion  error 
and  the  maximum  number  of  generations  that  the  GA  will  run.  For  the  ISAR  parameter 
search,  the  initial  population  is  a  matrix  with  the  first  element  of  each  parent  being  f  as 

determined  by  the  FFT.  (Note  that  from  this  point  on  when  search  parameters  are  referred 
to  with  a  single  subscribe,  the  argument  applies  equally  to  translational  and  rotational 
motion  compensation  search  parameters.  The  double  subscript  will  be  added  when  nec¬ 
essary  to  differentiate  between  the  two  sets  of  search  parameters.)  This  parameter  does 
not  change.  The  remaining  parameters,  {/2,/3,...} ,  are  initialized  as  uniform  random 

variables  between  ±N,  where  N  is  the  number  of  pulses  in  the  cross-range.  After  initiali¬ 
zation,  the  population  will  resemble  the  construct  shown  in  Table  1. 
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Degree 


P 
O 
P 
U 
L 
A 
T 
I 

O 
N 

S 
I 

z 

E 

V 

Table  1 .  The  Genetic  Algorithm  Population  -  2nd  subscript  indicated  Parent. 

b.  Evaluating  Fitness  of  Each  Member  of  Population 

In  any  search  algorithm,  the  ability  to  evaluate  the  suitability  of  a  set  of 
parameters  compared  to  other  sets  of  parameters  when  applied  to  the  cost  function  is  ab¬ 
solutely  crucial  in  order  to  ensure  that  the  best  parameter  set  is  chosen.  In  a  GA,  it  is 
critical,  since  the  fitness  ranking  will  determine  which  parents  will  be  used  to  form  the 
next  generation  and  which  of  the  newly  formed  children  will  be  inserted  into  the  new 
population.  Evaluating  fitness  is  the  key  component  in  determining  the  search  pattern  of 
the  solution  space.  Since  this  GA  uses  the  FFT  version  of  the  automatic  recognition  of 
focus  method  to  determine  parameter  fitness,  all  of  the  parents  or  children  are  evaluated 
against  the  cost  function  in  Equation  (28).  In  addition,  because  the  calculated  fitness  be¬ 
tween  superior  and  poor  parameters  are  small,  a  scaling  exponential  (chosen  from  em¬ 
pirical  observation  to  be  100)  was  added  to  the  fitness  calculation  in  order  to  increase  the 
separation  between  the  fitness  score  of  parameters.  Care  must  be  taken  when  choosing 
the  exponential.  Too  low  and  it  will  make  no  difference  to  the  efficiency  of  the  search 
but  set  it  too  high  and  the  GA  will  exit  part  of  the  solution  space  and  may  fail  to  find  the 
global  maximum  that  could  reside  there. 
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c.  Generating  Children 

Now  that  the  fitness  of  each  parent  has  been  computed,  the  next  step  is  to 
form  the  children  that  will  be  part  of  the  next  generation.  The  first  step  is  to  determine, 
based  on  the  fitness  of  each  parent,  the  likeliness  that  a  parent  will  be  mated  with  another 
parent  to  form  two  children.  This  is  done  using  the  roulette  wheel  method  as  introduced 
in  [9]  and  [10].  The  way  the  roulette  wheel  works  is  that  each  parent  is  assigned  a  por¬ 
tion  of  the  wheel  that  is  proportional  to  the  percentage  they  contribute  to  the  sum  of  the 
fitness  values  of  the  entire  population.  This  way  parents  who  contribute  significantly  to 
the  population  have  a  significant  chance  of  being  chosen  to  produce  children  while  even  a 
poor  contributor  retains  a  chance,  albeit  a  small  one,  to  contribute  to  the  population.  It  is 
critical  that  even  poor  parents  be  given  a  chance  to  reproduce  since  it  ensures  genetic  di¬ 
versity  within  the  population  and  leads  to  a  complete  search  of  the  solution  space.  The 
roulette  wheel,  shown  in  Figure  19,  was  built  in  MatLab  by  assigning  a  vector  a  cumula¬ 
tive  sum  of  the  fitness  values  (i.e.,  first  value  equal  3,  second  value  equal  to  13.5  and  the 
last  equal  to  30.02). 

Parent  6 
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Figure  19.  Illustration  of  the  Roulette  Wheel. 
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Once  the  roulette  wheel  has  been  constructed,  the  mating  pool  can  be  con¬ 
structed.  As  in  Figure  19,  if  our  population  is  of  size  6,  we  need  6  parents  to  form  the 
mating  pool.  Choosing  is  done  by  generating  a  pointer  which  is  computed  as  a  linear  dis¬ 
tributed  random  variable  between  0  and  ^  Fitness .  A  parent  is  chosen  if  it  is  the  first 

parent  with  a  roulette  wheel  value  greater  than  the  random  number  assigned  to  the 
pointer.  Figure  20  shows  a  pointer  which  possesses  a  value  of  19.  On  the  roulette  wheel, 
Parent  3  has  a  value  of  1 8  while  Parent  4  has  a  value  of  2 1 .  Since  1 9  >  1 8  and  1 9  <  2 1 , 
Parent  4  is  selected  for  the  mating  pool.  Parent  selection  for  the  mating  pool  continues  in 
this  manner  until  all  6  slots  in  the  mating  pool  are  filled. 
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Pointer  =  19  In  this  case,  parent  4  is  selected 
Figure  20.  Illustration  of  the  Mating  Pool  Selection  Method. 

Once  the  mating  pool  has  been  populated  it  is  broken  into  two  equal  arrays 
and  their  order  is  randomized.  The  matching  elements  of  the  two  breeding  pools  are  the 
mating  pairs.  Breeding  between  the  two  mating  pairs  are  determined  by  the  cross-over 
equations  taken  from  [6]: 

Cj  =  aPxx  +  (1  -  a)Pl  2  (34) 

and 

C2  =  (1  -  a)Pu  +  aPl2  (35) 
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where  Pl  l  and  Pl  2  are  the  first  parents  in  the  two  breeding  pools,  C1  and  C2  are  their 

children  and  a  is  a  uniformly  distributed  random  variable  between  0  and  1 .  So  for  each 
search  parameter,  the  child  is  getting  a  new  parameter  that  is  equivalent  to 
L  =  (a)  f\  +(1  -a)L  .  Also  remember  that  f.  does  not  change  as  it  is  de- 

U  z CHILD]  V  u  z PARENT  1  V  U  Z PARENT 2  ^  1  ^ 

termined  by  the  FFT. 

d.  Mutation  of  the  Children 

The  purpose  of  mutation  is  to  allow  the  search  to  be  global  in  nature  even 
as  the  parents  converge  towards  each  other  and  a  possible  global  maximum.  The  muta¬ 
tion  process  occurs  right  after  the  children  have  been  formed  and  is  illustrated  in  Figure 
21.  Mutations  are  determined  by  [6]: 

Cx=aCx+(\-a)Pmndom.  (36) 

In  Equation  (36),  C,  is  the  child  to  be  mutated,  Pmndom  is  a  parent  pulled  randomly  from 

the  solution  space  and  a ,  as  before,  is  a  uniform  random  variable  distributed  between  0 
and  1. 


Mutation 


Solution 

Space 


Figure  2 1 .  Illustration  of  the  Mutation  Process. 
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The  other  parameter  that  needs  to  be  addressed  is  the  mutation  rate.  Set¬ 
ting  the  rate  too  high  destroys  the  integrity  of  the  local  search  of  a  maximum  and  the  GA 
may  be  terminated  before  it  is  able  to  find  the  absolute  maximum.  However,  setting  it 
too  low  increases  the  probability  that  the  GA  will  converge  to  a  local  maximum.  Refer¬ 
ence  [10]  recommends  a  low  mutation  rate  of  0.1%  -  1%  since,  in  nature,  mutations  are 
relatively  rare.  Reference  [6]  uses  a  mutation  rate  of  30%  for  their  ISAR  GA.  Reference 
[11]  suggests  a  variable  rate,  one  which  starts  high  and  reduces  during  the  number  of  cy¬ 
cles.  Empirically,  it  was  determined  that  the  mutation  rate  from  [6],  set  at  a  constant  30% 
,  produced  the  best  results.  In  this  thesis’  GA,  this  mutation  rate  was  required  to  reliably 
break  the  GA  out  of  a  situation  where  it  would  repeatedly  converged  to  a  local  maximum 
that  created  a  poorly  focused  image. 

e.  Evaluating  the  Fitness  of  Each  Child  and  Reinsertion  of  Parents 

and  Children  into  New  Population 

Once  the  mutation  process  has  been  completed,  each  child  is  evaluated  for 
its  fitness  as  described  above.  The  next  stage  is  to  devise  a  way  to  reinsert  the  old  parents 
and  the  new  children  into  the  new  population  which  promotes  convergence  to  the  global 
maximum.  Reference  [10]  suggests  that  reinsertion  should  be  a  random  process  where 
only  a  certain  percentage  of  the  children  are  randomly  chosen  to  join  the  new  population 
and  the  remaining  empty  slots  are  randomly  filled  by  parents  of  the  old  population.  Care 
should  be  taken  when  applying  this  approach,  particularly  when  speed  of  convergence  is 
of  particular  importance.  What  occurs  in  this  approach  is  that  the  current  maximum  (i.e., 
the  best  value  found  so  far  in  the  search)  could  depart  the  population,  poor  performers 
will  enter  into  the  population  and  high  ranked  children  may  not  enter  the  population.  A 
compromise  to  this,  which  was  initially  implemented  in  this  thesis’  GA,  was  setting  the 
reinsertion  rate  of  the  children  to  an  80%  probability  of  insertion  into  the  new  population 
with  the  remaining  empty  slots  in  new  population  being  filled  by  the  highest  ranked  par¬ 
ents.  The  probability  of  convergence  graph,  as  a  function  of  size  of  population  and  the 
number  of  generations  evolved  using  this  reinsertion  method,  was  generated  and  is  dis¬ 
played  in  Figure  22.  (Discussion  on  how  the  convergence  graph  is  generated  will  be  dis¬ 
cussed  later  in  the  results  section.)  As  can  be  seen  from  this  graph,  the  probability  of 
convergence,  for  30  trials  at  each  point,  has  a  disappointing  maximum  value  of  80%  even 

with  a  large  number  of  generations  and  a  large  population  size.  To  improve  this,  it  is 
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tempting  to  choose  a  deterministic  approach  where  the  best  parents  and  children  are  cho¬ 
sen  to  form  the  new  population.  However,  while  building  the  GA  for  this  thesis,  imple¬ 
menting  reinsertion  in  this  manner  tended  to  lead  to  convergence  to  a  local  maximum  and 
not  the  global  maximum. 


Percentage  of  Time  Genetic  Algorithm  Converged  to  Focused  Image 
B727  Data  Set  -  30  Trials  Each  Point 


Figure  22.  Convergence  Graph  of  GA  with  Probability  of  Reinsertion  of  Children 

into  New  Population  Equal  to  80% 

In  order  to  improve  the  probability  of  convergence,  a  simple  method 
which  combines  deterministic  and  random  processes  was  chosen.  When  a  new  popula¬ 
tion  is  ready  to  be  formed,  50%  of  the  slots  are  filled  by  the  best  parents  of  the  old  popu¬ 
lation  and  the  best  children  of  the  new  population.  The  remaining  slots  are  filled  by  ran¬ 
domly  choosing  from  the  parents  and  children  not  yet  chosen  for  insertion  into  the  new 
population.  This  method  accomplishes  two  goals.  First,  it  ensures  that  the  best  values 
found  by  the  GA  will  continue  to  influence  the  evolution  of  the  GA  from  generation  to 
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generation  which  encourages  local  searches  around  these  good  parameter  sets.  Secondly, 
the  random  selection  of  the  second  half  of  the  new  population  ensures  sufficient  genetic 
diversity  in  the  population  which  encourages  global  searching  of  the  entire  solution 
space.  This  is  the  method  that  was  employed  in  this  thesis’  GA.  Another  convergence 
graph  was  generated  with  this  change  but  it,  and  its  subsequent  improvement,  are  dis¬ 
cussed  in  the  results  section. 

f.  Determine  if  the  Search  is  Complete 

There  are  a  couple  of  ways  to  determine  if  the  GA  search  is  complete. 

The  first  is  a  deterministic  way  where  the  user  indicates  the  maximum  number  of  genera¬ 
tions  that  the  search  algorithm  is  allowed  to  evolve.  After  the  last  generation,  the  best  pa¬ 
rameter  is  taken  and  used  for  either  translational  or  rotational  motion  compensation.  The 
second  method  is  to  measure  the  number  of  generations  that  the  best  value  found  has  re¬ 
mained  unchanged.  The  longer  this  value  remains  stagnant,  the  likelihood  of  finding  a 
better  value  becomes  less  and  less.  After  a  fixed  number  of  generations  that  this  value 
remains  stagnant,  the  GA  will  stop  and  use  this  value  for  motion  compensation.  For  this 
thesis,  the  user  indicates  the  number  of  generations  the  GA  will  run  since  this  allows  for 
the  generation  of  convergence  graphs,  which  in  turn  is  a  good  measure  of  the  GA’s  per¬ 
formance  and  allows  the  GA  to  be  compared  to  the  performance  of  the  PSO  algorithm. 

Finally,  this  procedure  must  be  completed  twice,  once  for  translational 
and  then  for  rotational  motion  compensation.  If  compared  to  the  exhaustive  search,  the 
number  of  cost  function  calculations  for  a  GA  can  be  determined  by: 


2. 


N  = 2- N  ■  N 

cos  t  _  function  _  GA  pop  Generations  * 


GA  ISAR  Results 


(37) 


To  determine  the  functionality  of  the  AJTF  GA,  the  B727  simulated  data  set  and 
several  different  2048  pulse  imaging  intervals  of  the  6  -  Point  Delta  Wing  data  set  were 
used. 


a.  B727  Simulated  Data  Set 

Shown  is  Figure  23  is  the  unfocused  B727  and  the  power  spectras  of  the 
two  range  bins  used  in  the  AJTF  GA.  It  is  easy  to  see  from  both  the  image  and  the  power 
spectras  that  there  is  time-varying  Doppler  in  the  simulated  radar  data  set.  Figure  24 
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shows  the  image  after  translational  motion  compensation  of  3rd-degree  motion.  Since 
most  of  the  error  is  rotational,  as  stated  in  [3],  the  image  is  still  not  focused  with  the  ex¬ 
ception  of  one  of  the  two  point  scatterers  in  range  bin  31,  which  was  changed  in  its  spec¬ 
trum  from  being  very  broad  to  very  narrow.  Figure  25  shows  the  final  focused  image. 
Notice  that  the  image,  after  3  -degree  rotational  motion  compensation,  is  well  focused 
and  that  the  scatterers  are  now  sharp  points  in  the  frequency  domain.  This  image  was 
formed  with  a  population  size  of  20  and  was  evolved  for  20  generations  for  a  total  run 
time  of  0.801  seconds  for  a  Ncost  function  GA  =  800.  As  given  in  Equation  (33),  an  exhaus¬ 
tive  search  would  have  taken  2048  calculations  of  the  cost  function. 
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Figure  23.  Unfocused  B727  and  the  Power  Spectras  of  Range  Bins  31  and  44. 
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Down  Range  -  Range  Bins 


Image  After  Translational  Motion  Compensation 
B727  Data  Set 
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Figure  24.  B727  and  the  Power  Spectras  of  Range  Bins  31  and  44  after  Transla¬ 

tional  Motion  Compensation. 


Final  Image  after  Rotational  Motion  Compensation 
B727  Data  Set 
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Figure  25.  Focused  B727  and  the  Power  Spectras  of  Range  Bins  31  and  44. 
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n  order  to  better  understand  the  relationship  between  the  population  size,  the  number  of 
generations  and  the  probability  of  convergence,  a  convergence  graph  was  calculated.  It 
was  built  by  varying  the  population  size,  increasing  the  number  of  generations  that  the 
GA  would  cycle  through  and,  after  50  tests  for  focus  had  been  made  for  each  point, 
measuring  the  percentage  of  times  the  unfocused  image  came  into  focus.  Focus  is  meas¬ 
ured  by  taking  the  2D  cross-correlation  of  the  focused  B727  image  with  the  resulting  im¬ 
age  generated  by  the  GA.  Cross-correlating  an  unfocused  image  with  the  focused  image 
will  create  a  smeared  2D  cross-correlation  plane  with  several  points  in  red  as  shown  in 
Figure  26.  Cross-correlating  a  focused  image  with  the  focused  image  will  create  a  2D 
cross-correlation  plane  with  one  distinct  peak  as  in  Figure  27.  This  can  easily  be  detected 
through  threshold  detection. 


Cross-Correlation  Function  of  an  Unfocused  and  Focused  B727  Image 
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Figure  26.  Cross-Correlation  between  a  Focused  and  Unfocused  B727  Image. 
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Cross-Correlation  Function  of  Two  Focused  B727  Images 
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Figure  27.  Correlation  between  Two  Focused  B727  Images. 


Figure  28  shows  the  convergence  graph  as  the  percentage  of  images  fo¬ 
cused  as  a  function  of  varying  the  population  size  and  the  number  of  generations  the  al¬ 
gorithm  was  allowed  to  evolve.  Figure  29  shows  the  average  run  time  of  each  of  these 
trials.  These  graphs  are  very  useful  in  determining  optimal  search  parameters  for  the  GA 
search.  For  example,  from  these  graphs  it  can  be  seen  that  an  optimal  choice  of  parame¬ 
ters,  for  the  B727  simulated  data  set,  would  be  a  GA  search  with  population  size  of  60 
which  is  evolved  for  10  generations  since  this  gives  an  84%  chance  of  convergence,  a  run 
time  of  only  0.96  seconds  and  only  1200  calculations  of  the  cost  function.  The  require¬ 
ment  to  calculate  the  cost  function  1200  times  is  only  58.5%  of  the  minimum  cost  func¬ 
tion  calculations  needed  by  the  exhaustive  search. 
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Percentage  of  Trials  the  Genetic  Algorithm  Converged  to  Focused  Image 
B727  Data  Set  -  50  Trials 


Figure  28.  Percentage  of  Images  Focused  Using  the  Genetic  Algorithm  as  a  Func¬ 
tion  of  Population  Size  and  Number  of  Generations  -  B727  Data  Set. 
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Average  Run  Time  of  the  Genetic  Algorithm 
B727  Data  Set  -  50  Trials 


Figure  29.  Average  Run  Time  of  the  Genetic  Algorithm  as  a  Function  of  Population 
Size  and  Number  of  Generations  -  B727  Data  Set. 

b.  6-  Point  Delta  Wing  Experimental  Data  Set 
The  6  -  Point  Delta  Wing  experimental  data  set,  with  2048  (211)  pulses  in 
the  cross-range,  represents  a  more  challenging  focusing  problem  as  it  does  not  always 
exhibit  ideal  motion  error  characteristics  like  the  B727  simulated  data  set.  Though  the 
assumptions  of  rigid  body,  two-dimensional  motion  still,  for  the  most  part,  hold  for  this 
data  set,  there  are  many  imaging  intervals  for  which  they  do  not  hold.  For  the  moment, 
the  discussion  of  results  will  be  limited  to  four  imaging  intervals  and  only  2nd-order  mo¬ 
tion  error  (i.e.,  a  1 -dimensional  solution  space  search  for  parameter  fn  and  f22 )  where 
the  GA  is  capable  of  finding  a  parameter  that  will  focus  the  image.  For  comparison,  the 
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number  of  cost  function  calculations  for  translational  and  rotational  motion  compensation 
using  the  exhaustive  search,  as  in  Equation  (32),  is  8196  calculations  of  the  cost  function. 

The  first  imaging  interval  is  imaging  interval  15  and  is  shown  in  Figure 
30.  This  imaging  interval  has  only  one  focused  point  scatterer  which  can  be  placed  in  its 
proper  cross-range  bin  by  the  fast  Fourier  transform.  The  remaining  scatterers  possess 
time-varying  Doppler  and  are  badly  blurred  in  the  cross-range.  The  GA  AJTF  algorithm 
with  a  population  of  size  30  can  find  the  search  parameter  to  focus  the  image  after  50 
generations  which  equates  to  3000  calculations  of  the  cost  function  or  37%  of  the  calcula¬ 
tions  when  compared  to  the  exhaustive  search.  Notice  that  the  power  spectrum  in  the 
unfocused  image  is  spread  in  Doppler  while  the  focused  image  has  a  power  spectrum 
with  sharp  points. 

The  second  imaging  interval  is  interval  19  and  is  shown  in  Figure  31. 
Again,  five  of  the  scatterers  are  out  of  focus  in  the  cross-range.  The  AJTF  GA  is  success¬ 
ful  in  finding  the  first  order  parameters,  fl2  and  f22 ,  and  focusing  the  image.  All  six 

point  scatterers  can  be  clearly  identified  in  the  image.  Figure  32  shows  imaging  interval 
23  whose  unfocused  image  has  significant  motion  error  as  only  one  comer  reflector  is  in 
its  proper  Doppler  bin.  The  focusing  algorithm  removes  the  error  and  focuses  the  image 
in  the  cross-range.  The  sixth  reflector  is  located  in  range  bin  20.  However,  as  it  is  shad¬ 
owed  by  a  reflector  on  the  leading  edge  of  the  6  -  Point  Delta  Wing,  in  range  bin  26,  it  is 
below  threshold  and  not  visible.  Again  both  the  imaging  intervals  were  formed  with 
3000  calculation  of  the  cost  function. 

The  final  imaging  interval  to  be  presented  is  imaging  interval  14  in  Figure 
33.  The  original  image  is  very  badly  blurred  in  the  cross-range.  The  AJTF  GA  is  able  to 
focus  most  of  the  unfocused  scatterers  extremely  well.  There  still  is  some  blurring,  par¬ 
ticularly  in  range  bin  25  but  also  a  small  amount  in  range  bins  26,  21  and  28.  The  inabil¬ 
ity  of  the  AJTF  algorithm  to  remove  all  blurring  is  most  likely  due  to  a  violation  of  one  of 
the  basic  assumptions  used  to  formulate  the  mathematical  model  of  the  AJTF,  most  likely 
the  presence  of  a  very  small  degree  of  3D  motion. 
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Figure  30.  Uncompensated  (top)  and  Motion  Compensated  6  -  Point  Delta  Wing 
ISAR  Image  -  Image  Interval  15  with  2048  Pulses  in  the  Cross-Range  -  GA. 
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Uncompensated  Image 
Delta  Wing  Data  Set  -  Imaging  Interval  19 
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Figure  3 1 .  Uncompensated  (top)  and  Motion  Compensated  6  -  Point  Delta  Wing 
ISAR  Image  -  Imaging  Interval  1 9  with  2048  Pulse  in  the  Cross-Range  -  GA. 
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Down  Range  -  Range  Bins 
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Final  Image  after  Rotational  Motion  Compensation 
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Figure  32.  Uncompensated  (top)  and  Motion  Compensated  Delta  Wing  ISAR  Image 
-  Imaging  Interval  23  with  2048  Pulses  in  the  Cross-Range  -  GA. 


52 


Uncompensated  Image 
i  Wing  Data  Set  -  Imaging  Interval  14 
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Final  Image  after  Rotational  Motion  Compensation 
Delta  Wing  Data  Set  -  Imaging  Interval  14 
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Figure  33.  Uncompensated  (top)  and  Motion  Compensated  Delta  Wing  ISAR  Image 
-  Imaging  Interval  14  with  2048  Pulses  in  the  Cross-Range  -  GA. 
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As  with  the  B727  data  set,  a  convergence  graph  for  the  Delta  Wing  data 
set  was  calculated  in  order  to  determine  the  best  combination  of  population  size,  number 
of  generations  in  the  run  time  and  the  resulting  probability  of  convergence.  Imaging  in¬ 
terval  15  was  used  as  the  test  data  set  to  be  focused  by  the  GA.  The  resulting  image  was 
then  compared  to  the  focused  version  of  this  imaging  interval  via  the  2D  cross-correlation 
function  following  the  methodology  as  explained  for  the  B727  data  set. 

The  convergence  graph  is  shown  in  Figure  34.  From  the  graph,  it  can  be 
seen  that  for  an  80%  probability  of  convergence,  a  population  size  of  50  and  a  run  of  50 
generations  is  necessary.  This  translates  into  5000  calculations  of  the  cost  function,  or 
61%  of  what  is  necessary  for  the  exhaustive  search.  It  is  shown  in  Figure  35  that  the  av¬ 
erage  run  time  for  a  GA  with  this  population  size  and  50  generational  evolutions  is  22 
seconds. 


Percentage  of  Trials  the  Genetic  Algorithm  Converged  to  Focused  Image 
Delta  Wing  Data  Set  -  2nd  Order  Motion  -  50  Trials 


Figure  34.  Percentage  of  Images  Focused  Using  the  Genetic  Algorithm  as  a  Func¬ 
tion  of  Population  Size  and  Number  of  Generations:  6  -  Point  Delta  Wing  Data 

Set. 
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Average  Run  Time  of  Genetic  Algorithm 
Data  Set-  2nd  Order  Motion  -  50  Trials 


Figure  35.  Average  Run  Time  of  the  Genetic  Algorithm  as  a  Function  of  Population 
Size  and  Number  of  Generations:  6  -  Point  Delta  Wing  Data  Set. 

C.  PARTICLE  SWARM  OPTIMIZATION  PARAMETER  SEARCH 

The  particle  swarm  optimization  (PSO)  algorithm  was  introduced  in  [7]  and  finds 
its  roots  in  the  fact  that  the  way  that  fish  school  and  bees  swarm  is  optimal  in  nature  and 
can  be  adapted  for  parameter  searches  to  optimally  solve  engineering  problems.  In  the 
case  of  the  ISAR  focusing  problem,  it  is  ideally  suited  to  finding  the  optimal  values  for 
the  search  parameters,  {/12,/B,...}  and  {/22,/23,--} ,  as  it  displays  the  same  search  char¬ 
acteristics  as  GA,  in  particular,  the  ability  to  rapidly  converge  towards  a  global  maxi¬ 
mum.  Figure  17  and  the  arguments  that  accompany  it  about  initial  rapid  improvement  to 
the  found  global  maximum  followed  by  smaller  improvements  after  many  cycles  applies 
equally  to  the  PSO  as  explained  for  the  GA.  The  driving  force  behind  the  particle  swarm 
optimization  is  that  each  particle  in  the  swarm  knows  the  fitness  of  its  current  location, 
the  best  location  it  has  ever  encountered,  the  Local  Best,  and  the  best  location  ever  en- 
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countered  by  any  particle  of  the  swarm,  the  Global  Best.  It  is  this  knowledge  that  forms 
the  search  pattern  of  the  solution  space  in  order  to  find  the  optimal  solution. 

1.  Sequence  of  the  Particle  Swarm  Optimization  Algorithm 
Like  the  GA,  the  PSO  follows  a  set  sequence  that  repeats  until  done.  Figure  36 
graphically  describes  each  of  the  steps  in  the  PSO  as  it  applies  to  ISAR  focusing.  The 
following  paragraphs  describe  each  step  in  detail. 


Figure  36.  The  PSO  Flowchart  (After  [7].) 
a.  Defining  the  Particles  of  the  Swarm 

The  swarm  in  a  PSO  is  identical  in  definition  to  the  population  in  a  GA,  as 
shown  in  Table  1 .  Each  particle  in  the  swarm  is  analogous  to  a  parent  and  holds  the 
search  parameters  {fvf2,fr..} ,  as  determined  by  the  degree  of  motion  error.  Each  parti¬ 
cle  is  initialized  as  a  uniform  random  variable  from  +  N  to  -N ,  where  N  is  the  number 
of  pulses  in  the  cross-range.  At  the  same  instance  of  defining  the  swarm,  each  particle  is 
designated  an  initial  velocity  that  is  uniformly  random  between  -2N  and  +2 N.  The 
maximum  velocity,  |  vmax  |,  that  any  particle  can  possess  is  2 N.  Figure  37  shows  a 

swarm  of  a  3rd-degree  parameter  search  after  initialization.  As  expected,  it  is  character¬ 
ized  by  the  uniform  spread  of  the  particles  through  the  expanse  of  the  solution  space. 
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Initial  Swarm  with  Particles 

Randomly  Scattered  in  a  3D  Solution  Space  -  4th-Degree  Motion  Error 


Figure  37.  Example  of  an  Initial  Swarm  in  a  Solution  Space  for  4th-Degree  Motion 

Error. 

b.  Assessing  Initial  Fitness  of  the  Swarm  and  Defining  the  Local 
Best  Vector  and  Global  Best  Particle 

Now  that  the  initial  swarm  has  been  defined,  each  particle  is  assessed  for 
its  fitness  as  defined  by  the  cost  function  in  Equation  (28)  in  a  process  that  is  identical  to 
that  described  in  the  GA  section  of  this  chapter.  In  fact,  the  MatLab  code  that  defines  the 
fitness  is  identical  for  the  GA  and  the  PSO.  Once  the  fitness  is  known,  the  best  particle  in 
the  swarm  becomes  the  Global  Best  which  holds  the  location  of  this  particle  and  its  fit¬ 
ness  score.  The  Local  Best  vector  is  initialized  with  the  current  particle  positions  and 
their  current  fitness  since  this  position  is  each  particle’s  first  position  and,  therefore,  by 
default  each  particle’s  best  position  they  have  ever  encountered. 

c.  Updating  the  Particle  Velocities  and  Positions 

Now  that  each  of  particles  has  a  fitness  score  assigned  to  it  and  its  velocity 
initialized,  the  particle  velocities  can  be  updated  based  on  the  current  location  of  each  in¬ 
dividual  particle,  its  Local  Best  particle  and  the  Global  Best  particle  of  the  swarm.  To 
do  this  Equation  (3)  in  [7]  is  used: 


57 


V 


K (v„  +  (px  •  rand  •  (/ioca/>„  -  /„ )  +  (p2  •  rant/  •  (/GZoW>n  -  /„ )) 


(38) 


where  vn  is  the  particle’s  current  velocity,  fn  is  the  particle’s  current  location  in  the  nth 
dimension  of  the  solution  space  and  fLocal  n  and  fGlobal  „  are  the  Local  and  Global  best  po¬ 
sitions,  respectively,  in  the  n,h  dimension.  The  PSO  parameters  <px ,  <p2  and  K  are  pa¬ 
rameters  that,  as  stated  by  [7],  were  determined  by  extensive  research  to  be  optimal  when 
set  to  2.8,  1.3  and  0.729,  respectively.  The  values  of  <px  and  <p2  controls  the  balance  be¬ 
tween  the  attraction  of  the  particle  to  the  local  and  global  bests.  With  <px>  cp2,  the  parti¬ 
cle  is  more  prone  to  be  attracted  towards  searching  the  area  surrounding  the  particle’s  lo¬ 
cal  best.  The  constriction  factor,  K ,  can  be  compared  to  inertia  and  it  effectively  assists 
in  the  deceleration  of  the  particle  as  convergence  occurs.  Rand  is  a  uniformly  distributed 
random  variable  between  0  and  1  that  serves  the  purpose  of  adding  randomness  to  this  at¬ 
traction  and  leads  to  a  more  complete  spanning  of  the  solution  space.  We  can  see  from 
Equation  (38)  that,  as  the  separation  between  the  particle’s  current  position  and  the 
swarm’s  Global  Best  and  its  Local  Best  expands,  the  greater  effect  they  will  have  on  the 
particle’s  velocity  and  subsequently  its  search  direction.  Also,  the  search  pattern  of  the 
entire  swarm  is  instantaneously  affected  every  time  the  Global  Best  changes.  Reference 
[7]  states  that  this  equation  is  analogous  to  a  bee’s  desire  to  keep  searching  for  food  in  the 
direction  that  it  is  going  but  being  instinctively  pulled  towards  the  best  location  it  has 
ever  found  (Local  Best)  and  the  best  location  ever  found  by  any  member  of  the  swarm 
(Global  Best). 

To  update  the  particle’s  positions,  it  is  simply  a  matter  of  applying  Equa¬ 
tion  (39)  to  each  of  the  search  parameters  [7]: 

fn=fn+vnt.  (39) 

Typically  t  is  taken  to  be  equal  to  1  so,  to  generate  the  next  cycle  of  the  swarm,  it  simply 
becomes  a  matter  of  adding  the  velocities,  vn ,  to  each  appropriate  search  dimension’s 

current  positions.  The  effect  that  the  velocity  update  has  on  the  particles  is  shown  in  Fig¬ 
ure  38. 
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Figure  38.  Graphical  Illustration  of  the  Velocity  and  Position  Update  Process.  Note 
that  gbest  =  Global  Best  and pbest  =  Local  Best  (From  [7].) 

When  discussing  velocity  and  position  updates,  it  is  necessary  to  discuss 
the  boundary  conditions  of  the  solution  space  since  particles  left  to  themselves  will  depart 
the  solution  space  where  meaningful  solutions  reside.  Since  we  know  from  previous 
discussion  that  our  search  parameters,  {/2,/3,...} ,  will  possess  values  between  ±V  ,  each 

dimension  of  the  solution  space  has  a  boundary  at  ±N  (where  N  is  again  the  number  of 
pulses  in  the  cross-range).  Reference  [7]  states  that  there  are  three  ways  to  deal  with  par¬ 
ticles  that  exceed  the  boundary  of  the  solution  space: 

i.  Absorbing  Walls:  When  a  particle  exceeds  a  boundary  in  a  dimen¬ 
sion  of  the  solution  space,  the  velocity  in  this  dimension  is  immediately  zeroed  and  its 
position  is  set  to  the  boundary.  The  particle  will  then  rejoin  the  swarm  when  its  velocity 
is  updated  during  the  next  cycle.  Extreme  care  must  be  taken  when  using  this  boundary 
condition  since  a  serious  problem  that  is  not  mentioned  in  [7]  may  occur.  If  the  situation 
exists  where  a  local  maximum  exists  at  the  boundary  and  if  this  position  is  simultane¬ 
ously  the  Local  Best  and  the  Global  Best  position,  the  velocity  will  remain  stuck  at  zero 
and  the  particle  will  stay  fixed  to  the  boundary.  A  condition  was  observed  where  all  the 
particles  in  a  swarm  became  stuck  to  the  boundary  and  the  swarm  did  not  converge  to  the 
true  global  maximum  and  the  focused  image  after  motion  compensation  was  very  poor. 

ii.  Invisible  Walls:  Here  the  idea  is  to  let  the  swarm  cycle  without  in¬ 
terruption  and  any  particle  that  strays  beyond  the  boundaries  of  the  solution  space  is  not 
evaluated  against  the  cost  function.  Instead  it  is  immediately  assigned  a  fitness  equal  to 


rA 


-  -  original  velocity 

. velocity  toward  gbest 

—  velocity  toward  pbest 
- resultant  velocity 


gbest 


pbest  2 

ik 


pbest  i 


59 


zero.  As  such,  this  particle  outside  the  boundaries  will  not  encounter  any  new  local  or 
global  maximums.  This  will  cause  the  velocity  update  equation,  as  in  Equation  (38),  to 
lead  the  particle  back  into  the  solution  space  without  interruption  to  the  particle’s  natural 
search  pattern. 

iii.  Rebounding  Walls:  When  a  particle  passes  a  boundary,  its  posi¬ 
tion  is  immediately  set  to  the  boundary  and  it  is  given  a  velocity  that  points  away  from 
the  boundary.  This  new  velocity  is  equal  to 

v ,  =±0.1-  rand  •  |vmax  _  (40) 

where  vmax  is  the  maximum  velocity  of  a  particle  in  the  n,h  direction,  rand  is  the  usual 

uniformly  distributed  random  variable  between  0  and  1  and  the  direction  is  set  such  that  it 
is  heading  away  from  the  boundary.  The  10%  scaling  factor  keeps  the  particle  near  the 
boundary,  since  that  is  where  it  was  headed,  while  still  allowing  the  particle  to  retain 
momentum  and  avoid  the  problem  of  sticking  to  the  boundary,  as  explained  above  in  the 
rebounding  walls  boundary  condition. 

Reference  [7]  goes  into  greater  detail  on  the  strengths  and  weaknesses  of 
each  boundary  condition.  However,  the  PSO  algorithm  written  for  this  thesis  uses  the  re¬ 
bounding  walls  option  since  it  keeps  all  the  particles  within  the  solution  space.  This  cre¬ 
ates  the  condition  where  each  particle  will  have  a  new  position  to  evaluate  against  the 
cost  function  for  each  cycle,  a  critical  factor  when  speed  of  convergence  is  important. 
Also,  there  is  no  danger  of  a  problem  with  the  boundary,  as  explained  for  the  absorbing 
wall  boundary  condition. 

Finally,  when  updating  the  velocity  of  the  particles,  it  is  necessary  to  en¬ 
sure  that  no  particle  exceeds  |vmax  |  in  either  the  positive  or  negative  direction.  Reference 

[7]  states  that  vmax  has  been  determined  to  be  optimal  when  set  to  equal  the  dynamic 
range  of  the  solution  space  in  a  dimension.  This  is  why,  as  previously  stated, 

|vmax  |  =  2  •  N.  When  a  particle  is  observed  exceeding  the  maximum  velocity,  it  is  simply 

set  equal  to  ±vmax ,  with  the  direction  of  velocity  preserved. 
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d.  Updating  Particle  Fitness  and  Reassigning  Local  Best  Vector 
and  Global  Best  Particle 

Now  that  the  new  positions  of  each  of  the  particles  have  been  defined, 
their  current  fitness  values  are  assessed  as  previous  discussed.  If  the  fitness  score  of  any 
of  the  particles  exceeds  that  of  the  Global  Best,  it  becomes  the  new  Global  Best.  In  addi¬ 
tion,  each  of  the  particles  current  fitness  is  assessed  against  their  individual  Local  Best. 

If  the  fitness  of  their  current  position  exceeds  their  Local  Best,  the  current  position  be¬ 
comes  the  Local  Best.  In  this  way,  as  new  Local  Bests  and  new  Global  Bests  are  found, 
the  search  pattern  of  each  particle  and  the  entire  swarm  is  updated  to  bias  their  search  to¬ 
wards  these  areas. 

e.  Termination  of  the  Search 

At  the  termination  of  the  search,  the  swarm  should  resemble  the  swarm 
shown  in  Figure  39  and  should  have  converged  to  the  solution  space’s  global  maximum. 
The  PSO  is  terminated  after  a  fixed  number  of  cycles  as  chosen  by  the  user,  which  is  the 
same  way  that  this  thesis’  GA  is  terminated.  As  with  the  GA,  terminating  the  search  in 
this  manner  adds  a  deterministic  quantity  to  a  random  search  which  allows  the  PSO  to  be 
evaluated  for  performance  and  compared  to  the  GA  in  terms  of  probability  of  conver¬ 
gence.  Similar  to  this  thesis’  GA,  on  completion  of  the  last  cycle,  the  Global  Best  parti¬ 
cle  is  taken  and  used  for  translational  or  rotational  motion  compensation  as  applicable. 
Also,  the  number  of  cost  function  calculations  of  the  PSO  algorithm  can  be  determined 
using  Equation  (37),  where  the  variable  names  have  been  altered  to  conform  to  the  con¬ 
ventions  of  a  PSO: 


N, 


cos  t  _  function  _  Swarm 


=  2  *  N  *N 

particles  cycles  * 


(41) 
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Final  Swarm  with  Particles 

Randomly  Scattered  in  a  3D  Solution  Space  -  4th-Degree  Motion  Error 


Figure  39.  Example  of  a  Final  Swarm  in  a  Solution  Space  for  4th-Degree  Motion 

Error. 

2.  PSO  ISAR  Results 

To  determine  the  functionality  of  the  AJTF  PSO  algorithm,  the  B727  simulated 
data  set  and  the  four  2048  pulse  imaging  intervals  of  the  6  -  Point  Delta  Wing  experi¬ 
mental  data  set,  were  used.  These  four  imaging  intervals  were  the  same  intervals  that 
were  used  in  the  GA  ISAR  results. 

a.  B727  Simulated  Data  Set 

A  focused  B727  image  is  found  in  Figure  40  with  the  spectra  of  the  range 
bins  used  for  focusing  shown  on  the  right  of  the  image.  It  was  focused  by  the  PSO  algo¬ 
rithm  in  0.681  seconds  after  20  cycles  and  with  a  swarm  size  of  20,  or  800  calculations  of 
the  cost  function.  As  with  the  GA,  it  is  much  more  important  how  this  algorithm  behaves 
in  a  consecutive  series  of  runs  rather  than  just  in  one  run  instance.  Using  the  PSO  algo¬ 
rithm,  a  convergence  and  run  time  graph  was  constructed  to  verify  algorithm  perform¬ 
ance  when  tasked  to  focus  this  simulated  data  set.  The  convergence  graph  is  shown  in 
Figure  41  while  the  average  run  time  graph  is  shown  in  Figure  42.  These  two  graphs 
show  that  an  optimal  selection  of  parameters  could  be  10  particles  and  50  cycles  as  this 
produces  an  88%  chance  of  convergence,  a  run  time  that  is  below  1  second  and  requires 
49%  of  the  required  cost  function  calculations  of  an  exhaustive  search.  If  probability  of 
convergence  was  of  critical  importance,  the  parameters  of  26  particles  and  30  cycles 
could  be  chosen.  These  parameters  yield  a  94%  chance  of  convergence  at  a  cost  of  a  run 
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time  of  1.2  seconds  and  76.1%  of  the  cost  function  calculations  as  compared  to  the  ex¬ 
haustive  search. 


Final  Image  after  Rotational  Motion  Compensation 
B727  Data  Set 
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Figure  40.  B727  Data  Set  Focused  by  the  PSO  Algorithm. 


Percentage  of  Trials  the  Particle  Swarm  Algorithm  Converged  to  Focused  Image 
B727  Data  Set  -  50  Trials 
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Figure  4 1 .  Percentage  of  Images  Focused  Using  the  Particle  Swarm  Optimization 
Algorithm  as  a  Function  of  Swarm  Size  and  Number  of  Cycles  -  B727  Data  Set. 
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Average  Run  Time  of  the  Particle  Swarm  Algorithm 
B727  Data  Set-  50 Trials 


Figure  42.  Average  Run  Time  of  the  Particle  Swarm  Optimization  Algorithm  as  a 
Function  of  Swarm  Size  and  Number  of  Cycles. 

b.  6-  Point  Delta  Wing  Experimental  Data  Set 
The  Particle  Swarm  Algorithm  was  tested  against  four  2048-pulse  imaging 
intervals  of  the  6  -  Point  Delta  Wing  experimental  data  set  in  order  to  ensure  that  it  could 
focus  real  radar  data.  These  four  imaging  intervals,  15,  19,  23  and  14,  are  the  same  inter¬ 
vals  that  were  used  in  the  GA  section  of  this  chapter  and  are  shown  in  Figures  43  -  46. 
The  focused  images,  which  were  generated  by  the  PSO  algorithm,  are  exactly  the  same  as 
those  generated  by  the  GA.  The  images  were  typically  formed  with  1200  calculations  of 
the  cost  function  over  a  run  time  of  6  seconds.  Again,  to  better  understand  the  relation¬ 
ship  between  swarm  size,  number  of  particles  and  their  effect  on  the  probability  of  con¬ 
vergence  and  run  time,  a  convergence  graph  and  a  run  time  graph  were  constructed  and 
are  displayed  in  Figures  47  and  48. 
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Final  Image  after  Rotational  Motion  Compensation 
Delta  Wing  Data  Set  -  Imaging  Interval  15 
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Figure  43. 


Motion  Compensated  6  -  Point  Delta  Wing  ISAR  Image  -  Image  Inter¬ 
val  15  with  2048  Pulses  in  the  Cross-Range  -  PSO. 


Final  Image  after  Rotational  Motion  Compensation 

Delta  Wing  Data  Set- Imaging  Interval  19  Power  Spectrum 
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Figure  44.  Motion  Compensated  6  -  Point  Delta  Wing  ISAR  Image  -  Image  Inter¬ 
val  19  with  2048  Pulses  in  the  Cross-Range  -  PSO. 
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Figure  45.  Motion  Compensated  6  -  Point  Delta  Wing  ISAR  Image  -  Image  Inter 
val  23  with  2048  Pulses  in  the  Cross-Range  -  PSO. 
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Figure  46.  Motion  Compensated  6  -  Point  Delta  Wing  ISAR  Image  -  Image  Inter 
val  14  with  2048  Pulses  in  the  Cross-Range  -  PSO. 
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From  Figures  47  and  48,  it  can  be  seen  that  the  PSO  algorithm  performs 
very  well  against  the  6  -  Point  Delta  Wing  experimental  data  set  and  for  an  84%  prob¬ 
ability  of  convergence,  it  requires  an  average  run  time  of  only  14.3  seconds  and  3000  cal¬ 
culations  of  the  cost  function.  This  represents  37%  of  the  cost  function  calculations  re¬ 
quired  by  the  exhaustive  search.  Also,  Figures  47  and  48  shows  that  the  PSO  is  capable 
of  performing  extremely  well,  even  with  a  small  number  of  particles  in  the  swarm.  For 
example,  using  20  particles  and  running  the  algorithm  for  50  cycles  leads  to  a  total  of 
2000  calculations  of  the  cost  function  (24%  of  the  exhaustive  search),  a  9-second  run 
time  and  a  98%  chance  of  convergence  to  focus.  This  shows  that  this  algorithm  pos¬ 
sesses  both  speed  and  reliability  even  against  a  real  radar  image  with  a  large  number  of 
pulses  in  the  cross-range. 


Percentage  of  Trials  the  Particle  Swarm  Algorithm  Converged  to  Focused  Image 
Delta  Wing  Data  Set  -  2nd  Order  Motion  -  50  Trials 


Figure  47.  Percentage  of  Images  Focused  Using  the  Particle  Swarm  Algorithm  as  a 
Function  of  Swarm  Size  and  Number  of  Cycles  -  Delta  Wing  Data  Set. 
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Average  Run  Time  of  Particle  Swarm  Algorithm 
Delta  Wing  Data  Set  -  2nd  Order  Motion  -  50  Trials 


Figure  48.  Average  Run  Time  of  the  Particle  Swarm  Optimization  Algorithm  as  a 
Function  of  Swarm  Size  and  Number  of  Cycles  -  Delta  Wing  Data  Set. 

D.  COMPARISON  BETWEEN  GA  AND  PSO  ISAR  FOCUSING 
ALGORITHMS 

Comparisons  between  the  GA  and  the  PSO  algorithms  are  made  in  [7]  and  [11] 
and  they  typically  find  that  both  evolutionary  search  techniques  produce  similar  results. 

In  particular,  [11]  states  that  during  some  optimization  trials,  which  in  this  case  was  op¬ 
timizing  the  complex  weights  of  an  antenna  array,  the  GA  outperforms  the  PSO  while  in 
others  the  PSO  outperforms  the  GA.  In  the  case  of  the  ISAR  focusing  problem,  for  imag¬ 
ing  intervals  where  the  AJTF  algorithm  was  a  valid  mathematical  model,  the  resulting  fo¬ 
cused  images,  produced  by  the  GA  and  PSO  algorithms,  were  of  equal  quality  as  both 
search  techniques  were  capable  of  finding  the  optimal  search  parameters.  When  a  poorly 
focused  image  was  produced  by  the  GA  or  PSO,  it  was  from  an  imaging  interval  of  the  6 
-  Point  Delta  Wing  experimental  data  set  that  possessed  qualities  that  were  in  severe  vio¬ 
lation  of  the  assumptions  of  the  AJTF  motion  compensation  model.  These  imaging  inter¬ 
vals  were  of  equivalent  poor  quality  regardless  of  whether  the  GA  or  the  PSO  was  used 
for  the  focusing  parameter  search.  Also,  by  examining  the  average  run  time  graphs  of 
both  the  GA  and  the  PSO,  it  can  be  concluded  that,  for  equivalent  calculations  of  the  cost 
function,  the  two  algorithms  are  equally  fast  at  processing  through  the  user-defined  num¬ 
ber  of  generations  or  cycles.  However,  if  the  convergence  graphs  for  the  B727  simulated 
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data  set  in  Figures  28  and  41  and  the  convergence  graphs  for  the  6  -  Point  Delta  Wing 
experimental  data  set  in  Figures  34  and  47  are  examined,  it  is  shown  that  there  is  a  sig¬ 
nificant  difference  in  performance  between  the  two  search  algorithms.  Figure  49  com¬ 
bines  part  of  the  GA  and  the  PSO  results  from  their  respective  convergence  graphs  with 
the  lines  with  triangles  denoting  GA  results  and  lines  with  circles  denoting  PSO  results. 
What  can  be  seen  from  the  graphs  in  Figure  49  is  that  the  PSO  written  for  this  thesis  con¬ 
sistently  outperforms  the  GA  by  reliably  converging  to  a  focused  image  even  when  the 
PSO  has  less  particles  in  the  swarm  than  the  GA  has  parents  in  its  population.  For  the 
B727  data  set,  a  PSO  with  10  particles  outperforms  a  GA  with  a  population  size  of  10  and 
16  and  performs  almost  as  well  as  a  GA  with  a  population  size  of  26.  In  the  6  -  Point 
Delta  Wing  experimental  data  set,  a  PSO  with  20  particles  outperforms  a  GA  with  a 
population  size  of  20,  30  and  50.  Also,  from  a  programming  point  of  view,  a  PSO  is 
much  easier  to  code  and  to  troubleshoot  and  its  search  pattern  is  judged  merely  by  Equa¬ 
tion  (38)  and,  although  some  of  the  constants  given  by  [7]  can  be  changed,  there  was  no 
requirement  to  implement  any  parameter  changes  for  the  PSO  in  this  thesis.  For  a  GA, 
setting  a  correct  mutation  rate,  reinsertion  rate  and  selecting  parents  correctly  for  breed¬ 
ing  are  all  of  critical  importance  for  reliable  convergence  and  can  take  considerable  effort 
to  set  at  a  correct  level  that  balances  global  and  local  searches.  For  these  reasons,  the 
PSO  is  determined  to  be  better  suited  for  ISAR  focusing  and  only  the  PSO  will  be  used 
for  the  subsequent  sections  of  this  thesis. 
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Percentage  of  Trials  the  Particle  Swarm  Algorithm  and 
Genetic  Algorithm  Converged  to  Focused  Image  -  B727  Data  Set-  50  Trials 
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Percentage  of  Trials  the  Particle  Swarm  Algorithm  and  Genetic  Algorithm 
Converged  to  Focused  Image  -  Delta  Wing  Data  Set  -  2nd  Order  Motion  -  50  Trials 


Figure  49.  Combined  GA  (top)  and  PSO  Convergence  Graphs  for  Selected  Popula¬ 
tion  /  Swarm  Sizes. 
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E.  HIGHER  ORDER  MOTION  COMPENSATION 

Higher  order  motion  occurs  when  a  target’s  motion  is  extremely  chaotic.  This 
creates  the  situation  where  searching  for  the  2nd-order  search  parameters,  fn  and  f22,  is 
insufficient  to  focus  the  ISAR  image.  This  has  been  partially  addressed  already  in  the 
B727  simulated  data  set  where  it  was  required  to  find  parameters  {/12,/13}  and 

{/22,/23}to  focus  the  image.  However,  with  only  256  pulses  in  the  cross-range,  which 

translates  to  an  ideal  solution  space  of  262,144  possible  combinations,  the  PSO  is  able  to 
quickly  traverse  the  solution  space  and  find  the  focusing  parameters.  The  problem  with 
higher-order  motion  compensation  occurs  when  there  are  a  large  number  of  pulses  in  the 
cross-range.  For  example,  for  3rd-order  motion  compensation,  the  solution  space  of  the 
2048  pulse  6  -  Point  Delta  Wing  experimental  data  set  increases  from  4096  possibilities 
to  16.7  x  106  possibilities.  If  4th-order  motion  compensation  is  required,  the  solution 
space  expands  to  68.7  x  109  possibilities.  Also,  do  not  forget,  that  these  solution  spaces 
must  always  be  traversed  twice,  once  to  extract  translational  and  the  other  to  extract  rota¬ 
tional  motion  compensation  parameters.  This  creates  an  enormous  search  space  that  can¬ 
not  be  quickly  traversed  by  even  the  most  efficient  PSO. 

Imaging  interval  13  of  the  6  -  Point  Delta  Wing  experimental  data  set  is  given  as 
an  example  and  its  unfocused  image  is  shown  in  Figure  50.  This  image  possesses  4th- 
order  motion  error  and  cannot  be  focused  with  2nd-order  motion  compensation.  The  exis¬ 
tence  of  higher-order  motion  can  be  seen  by  examining  a  range  bin  in  the  time-frequency 
plane.  As  explained  in  [3],  scatterers  with  higher-order  motion  will  be  curved  as  com¬ 
pared  to  2nd-order  motion  error  whose  scatterers  are  tilted  in  the  time-frequency  plane. 
The  time-frequency  transforms  of  the  scatterers  in  the  unfocused  range  bins  25  and  19  are 
shown  on  the  left  of  Figure  5 1 .  The  fact  that  they  are  curved  is  a  strong  indication  of  the 
presence  of  higher-order  motion  error.  If  the  standard  PSO  AJTF  focusing  algorithm,  us¬ 
ing  4th-order  motion  compensation,  is  applied  to  the  image,  an  extremely  well-focused 
image  is  produced  as  seen  in  Figure  52.  Also,  the  right  side  of  Figure  51  shows  that  the 
time-frequency  transforms  of  the  scatterers  from  the  focused  image  have  been  straight¬ 
ened  since  the  time-varying  Doppler  has  been  removed  from  the  signal.  Unfortunately, 
this  requires  approximately  1  x  106  calculations  of  the  cost  function.  This  limits  the  PSO 


71 


usefulness  since  this  requires  an  excessive  run  time  of  1  hour  and  10  minutes.  The  way 
around  this  is  to  adapt  the  AJTF  PSO  algorithm  to  perform  the  search  in  the  same  pattern 
as  the  exhaustive  search,  which  was  introduced  earlier  in  this  Chapter.  In  this  variation 
and,  like  in  the  exhaustive  search,  the  PSO  is  applied  to  searching  for  f2 .  Once  the  op¬ 
timal  fit  is  found,  /3  is  searched  for  with  the  PSO.  Once  the  optimal  f3  is  found,  as  long 
as  the  addition  of  the  /3  parameter  adds  to  the  fitness  value  as  calculated  by  Equation 
(28),  it  is  added  to  the  basis  function.  If  the  fitness  decreases,  /3  is  set  equal  to  zero. 

This  process  is  continued  for  as  many  degrees  of  motion  error  as  required  and  is  com¬ 
pleted  for  both  translational  and  rotational  motion  compensation.  The  resulting  focused 
image,  for  2nd,  3rd  and  4th-order  motion  compensation  is  shown  in  Figure  53.  The  4th- 
order  motion  compensated  image  is  shown  on  the  left  and  though  it  is  not  as  good  as  the 
image  in  Figure  52,  it  can  be  achieved  in  20  seconds  with  only  4500  calculations  of  the 
cost  function.  Therefore,  a  sub-optimal,  but  usable,  image  can  be  achieved  from  an  im¬ 
age  blurred  by  higher-order  motion  using  this  focusing  method  which  is  exceptionally 
fast  even  with  a  large  number  of  pulses  in  the  cross-range. 
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Figure  50.  Unfocused  Imaging  Interval  13. 
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Figure  5 1 .  Time-Frequency  Transforms  of  Range  Bins  25  and  19  -  Left:  Unfocused 

Right:  Focused. 
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Figure  52.  Focused  Imaging  Interval  13. 
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Figure  53.  Image  Interval  13  with  2nd,  3rd  and  4th-Degree  Motion  Compensation, 

Reduced  Solution  Space. 


F.  CHAPTER  SUMMARY 

In  this  chapter,  we  have  examined  the  construct  of  the  solution  space  in  the  ISAR 
focusing  problem  and  saw  that  it  possessed  many  local  maxima  but  only  one  global  maxi¬ 
mum.  This  characteristic  makes  the  solution  space  ill-suited  for  conventional  searches 
but  well  suited  for  the  GA  and  the  PSO  evolutionary  searches.  It  was  also  shown  that  the 
GA  and  the  PSO  provided  considerable  computational  time  saving  over  the  exhaustive 
search  which  is  the  method  normally  used  for  parameter  extraction  in  the  AJTF  algo¬ 
rithm.  Also,  it  was  shown  that  the  PSO  performed  considerably  better  than  the  GA  in  the 
ISAR  focusing  task.  Finally,  higher-order  motion  compensation,  using  the  PSO  AJFT, 
was  demonstrated  on  an  imaging  interval  of  the  6  -  Point  Delta  Wing  experimental  data 
set.  In  the  next  chapter,  we  will  examine  what  occurs  when  3D  motion  enters  the  radar 
signal  and  how  this  can  be  accommodated  since  this  is  a  violation  of  a  fundamental  as¬ 
sumption  of  the  AJTF. 
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V.  3D  MOTION  DETECTION 


In  Chapter  IV,  a  GA  and  PSO  algorithm  were  successfully  adapted  to  the  optimi¬ 
zation  of  the  AJTF  algorithm  which  provided  a  solution  to  the  computational  intensity 
required  for  successful  ISAR  focusing  that  was  one  of  the  algorithms  major  weaknesses. 
A  number  of  imaging  intervals,  which  are  now  well  focused  in  the  cross-range,  were  pre¬ 
sented.  These  imaging  intervals,  however,  follow  the  mathematical  model  of  the  AJTF. 
The  next  step  in  this  thesis  is  to  examine  an  effective  method  for  dealing  with  imaging  in¬ 
tervals  which  violate  the  AJTF’s  mathematical  model. 

A.  INTRODUCTION  TO  THE  3D  MOTION  MODEL  AND  THE 

DETECTION  METHOD 

One  of  the  basic  assumptions  of  the  adaptive  joint  time-frequency  motion  com¬ 
pensation  model  is  that  the  rotational  motion  used  to  form  the  ISAR  image  is  confined  to 
a  2D  plane  during  the  coherent  processing  interval,  TCPI .  When  motion  in  the  roll  plane 

of  the  aircraft  occurs  during  the  coherent  processing  interval,  as  shown  in  Figure  54,  3D 
motion  is  present  in  the  ISAR  return  signal  and  the  motion  error  cannot  be  compensated 
by  the  motion  compensation  model  in  its  current  form.  Due  to  the  difficulty  in  compen¬ 
sating  for  3D  motion  error,  Chen,  Li  and  Ling,  in  [3]  and  [4],  developed  the  linearity  of 
phases  method  to  detect  the  presence  of  3D  motion.  Thayaparan  further  expanded  on  the 
3D  detection  method  in  [12]  and  it  is  their  work  that  was  adapted  to  the  AJTF  PSO  mo¬ 
tion  compensation  algorithm  for  the  purposes  of  3D  motion  detection  in  this  thesis.  The 
goal  of  3D  motion  detection  is  to  determine  which  imaging  intervals  of  an  ISAR  data  set 
contains  a  low  degree  of  3D  motion  and  can  be  focused  using  the  AJTF  algorithm.  Imag¬ 
ing  intervals  found  to  possess  a  high  degree  of  3D  motion  are  simply  discarded  and  not 
presented  to  the  ISAR  operator. 
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Figure  54.  Engagement  of  a  Target  with  3D  Motion  (From  [3].) 

Now  that  there  is  3D  motion  in  the  ISAR  signal,  Equation  (1)  can  be  written  in  a 
more  general  form  of  [3]: 

s(x,t)  =  J^Ak  exp  |  -j^±[R(t)  +  xk  +  yk6(t)  +  ^(f)]|  (42) 

where  <j>(t)  describes  the  independent  angular  parameter  in  the  roll  plane  and  (xk,yk,zk ) 

are  the  3D  coordinate  position  of  the  kth  scatterer.  Chen  explains  in  [3]  that  Equation  (42) 
will  reduce  to  the  2D  model  in  two  cases.  The  first  case  is  when  zk  is  small  and  the 
phase  term  can  be  neglected.  The  second  case  is  when  there  exists  a  linear  relationship 
between  0{t)  and  tf>(t )  such  that  <f>(t)  =  a0(t)  [3]. 

The  first  step  towards  detecting  3D  motion  is  to  remove  translation  motion  error 
from  the  signal  by  finding  the  appropriate  basis  function,  h  (t) ,  as  previously  described 

and  multiplying  the  signal  by  its  conjugate.  With  the  translational  motion  error  removed, 
Equation  (42)  reduces  to  [4]: 
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(43) 


^  (x,t)T  =  J^Ak  exp|-y^^[Ax/(  +  Ayk0(t)  +  Az^(f)]  j 

=  YjA  exp{-y^^-[Axt  +  (Ay*  +  aAzk)0 (*)]' }• 

k=X  V  c  J 

Now  that  the  translational  motion  error  has  been  compensated,  it  becomes  necessary  to 
find  the  rotational  motion  compensation  search  parameters,  {/2i,/22,/23,...} ,  that  are  ap¬ 
propriate  to  formulate  the  phase  function,  0R  (?) ,  as  in  Equation  (20),  for  several  differ¬ 
ent  point  scatterers  in  the  target.  For  the  purposes  of  this  thesis,  3  point  scatterers  from 
the  6  -  point  Delta  Wing  data  set  are  used  and  the  AJTF  PSO  that  was  developed  above  is 
adapted  to  extract  the  search  parameters. 

Once  all  three  phase  functions  have  been  extracted  we  have: 

0)  =  fuZ  +  ^fn/2  +'- 

&Rl(t)  =  f2lj  +  —  fll/  +^j/232^3  +•••  (44) 

®«3  (0  =  fixf +Y\f22J2  +  + 

For  signals  that  do  not  possess  3D  motion,  any  of  the  above  phase  functions,  ©fil  it) , 

@R2  (t)  or  @R3  (t) ,  would  be  suitable  for  rotational  motion  compensation  which  would 

create  an  ISAR  image  that  is  well  focused  in  the  cross-range.  For  3D  motion  detection, 
Thayaparan’s  work  in  [12]  is  followed  and  an  average  phase  function  is  formed: 


0  (/)  = 

avg  V  L  / 


^  flX,  +  fl\2  +  fix 3  N 

i 

t  - i — 

^  f 22,  +fl22  +  f 22, ^ 

2  1 

t  - | - 

^  f  23,  +  f 232  +  f 23,^ 

V  3  J 

1  2! 

l  3  J 

1  3! ' 

v  3  J 

t  +... 


(45) 


The  purpose  to  formulating  0m,g  (t)  is  that  as  more  and  more  search  parameters  are  ex¬ 
tracted  from  more  and  more  point  scatterers  suitable  for  rotational  motion  compensation, 
@avg  (t) — m>maches  >  @jdeal  (t)  where  &ldeal  (t)  is  the  best  phase  function  possible  that  de¬ 
scribes  the  rotational  motion  error  in  the  whole  signal  and,  after  rotational  motion  com¬ 
pensation,  produces  the  best  image.  In  this  discussion,  the  number  of  point  scatterers  is 
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limited  to  three  since  the  maximum  number  of  scatterers  in  the  6  -  Point  Delta  Wing  data 
set  is  six  and  not  all  scatterers  are  suitable  for  rotational  motion  compensation. 


Now,  if  the  imaging  interval  possesses  a  low  degree  of  non-linearity,  the  plot  of 
&Rl,&R2  or  ®R3  as  a  function  of  ®ideal  would  be  a  straight  line  as  in  Figure  55.  There¬ 
fore,  the  first  step  in  determining  the  degree  of  non-linearity  is  to  find  the  straight  line  of 
best  fit,  in  a  least-squares  sense  using  the  MatLab  polyjit  function,  between  0ffl,0S2  or 

0A3  and  ®ideal .  The  second  step  is  to  plot  0AM  ,0/(l,  or  0A3  as  a  function  of  ®ideul  and 
measure  the  percentage  that  this  plot  deviates  from  the  line  of  least-squares  fit  using: 


N 


t=  1 


K,(e,A„(0)-LSF{e„(e, ,„,(;))}! 

LSF{0„,(0, „„,(())} 


(46) 


where  NF] \  is  the  degree  of  non-linearity  between  0AI  (t)  and  ®idea,  (t) ,  0ffl  (0i3fea/  (t))  is 
the  plot  of  0A1  {t)  as  a  function  of  ®jdeal  it) ,  LSF  {0AI  (®ideal  (>))}  is  the  straight  line 
least-squares  fit  between  0AI  (/)  and  ®ideal  (t)  and  N  is  the  number  of  pulses  in  the  cross¬ 
range.  This  procedure  is  repeated  to  generate  NF2  and  NF3  for  the  other  two  phase  func¬ 
tions  and  the  imaging  interval  is  assigned  its  degree  of  non-linearity  by  simply  averaging 
the  non-linearity  functions  [12]: 


NF, 


NF} ;  +  nf2  +  nf3 


Interval 


(47) 


Plot  of  ©m  as  a  Function  of  ©. .  . 
K1  ideal 


Figure  55.  Plot  of  0,,,  as  a  Function  of  ®ideal  when  no  3D  motion  is  present  (After 

[12].) 
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Figure  56  shows  the  Non-Linearity  of  Phases  plot  for  an  imaging  interval  with  a 
high  degree  of  non-linearity.  Note  that  the  plot  of  &Rl  ( Qi(leal  (t  ))  deviates  significantly 

from  the  line  of  least-squares  fit.  Figure  57  shows  the  Non-Linearity  of  Phases  plot  for 
an  imaging  interval  with  a  very  low  degree  of  non-linearity.  Note  that  the  plot  of 
©fll  (@lVfea/(0)  in  this  case  is  plotted  exactly  on  top  of  the  line  of  least- squares  fit. 


An  Imaging  Interval  that  Possesses  a  High  Degree  of 
Non  -  Linearity  of  Phases 


Figure  56.  Imaging  Interval  Possessing  a  High  Degree  of  Non-Linearity  (After 

[12].) 


An  Imaging  Interval  that  Possesses  a  Low  Degree  of 
Non  -  Linearity  of  Phases 


Figure  57.  Imaging  Interval  Possessing  a  Low  Degree  of  Non-Linearity  (After  [12].) 
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B.  3D  MOTION  DETECTION  RESULTS 

The  3D  motion  detection  algorithm  was  tested  against  the  6  -  Point  Delta  Wing 
data  set  with  2048  and  1024  pulses  in  the  cross-range  and  2nd  and  4th-degree  motion  com¬ 
pensation.  The  primary  goal  of  this  section  is  to  demonstrate  the  capability  of  Thaya- 
paran’s  3D  motion  detection  algorithm  to  differentiate  between  imaging  intervals  pos¬ 
sessing  a  high  degree  of  non-linearity  and  those  that  possess  a  low  degree  of  non-linearity 
after  the  algorithm  has  been  adapted  to  the  AJTF  PSO  algorithm,.  Secondly,  the  effect  of 
choosing  different  sizes  of  imaging  intervals  and  different  degrees  of  motion  compensa¬ 
tion  in  detecting  3D  motion  will  be  examined. 

1.  Imaging  Intervals  with  2048  Pulses  and  2nd-Degree  Motion  Compen¬ 
sation 

If  the  3D  motion  detection  algorithm  is  applied  to  the  6  -  Point  Delta  Wing  ex¬ 
perimental  data  set  with  its  imaging  intervals  divided  into  2048  pulses  and  2nd-degree 
motion  compensation  is  applied,  the  non-linearity  of  phases  graph  is  shown  in  Figure  58. 

Degree  of  Non-Linearity  in  Each  Imaging  Interval 
6  -  Point  Delta  Wing  Data  Set  with  2048  Pulses  and  Degree  of  Motion  Compensation  -  2 


Figure  58.  Degree  of  Non-Linearity  of  Imaging  Intervals  of  2048  Pulses  and 

2nd-Degree  Motion  Compensation. 
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Figure  58  indicates  that  imaging  intervals  with  a  low  degree  of  non-linearity  of 
phases  are  imaging  intervals  23,  15,  and  18.  If  these  images  are  focused  using  the  AJTF 
PSO  designed  for  this  thesis,  the  resulting  well  focused  images  are  shown  in  Figures  59  - 
61: 


Figure  59. 


Imaging  Interval  23  -  Well  Focused  with  Low  Degree  of  Non-Linearity 

of  Phases. 
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Figure  60.  Imaging  Interval  15  -  Well  Focused  with  Low  Degree  of  Non-Linearity 

of  Phases. 


Figure  61 .  Imaging  Interval  18  -  Well  Focused  with  Low  Degree  of  Non-Linearity 

of  Phases. 
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Figure  58  also  indicates  that  imaging  intervals  11,12  and  26  possess  a  high  de¬ 
gree  of  non-linearity  of  phases.  If  these  images  are  focused  using  the  AJTF  PSO  de¬ 
signed  for  this  thesis,  the  following  poorly  focused  images,  shown  in  Figures  62  -  64  re¬ 
sult: 


®ideal 


©..  . 
ideal 


Figure  62.  Imaging  Interval  1 1  -  Poorly  Focused  with  High  Degree  of  Non- 

Linearity  of  Phases. 


Imaging  Interval  12  High  Non-Linearity 
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Figure  63.  Imaging  Interval  12  -  Poorly  Focused  with  High  Degree  of  Non- 

Linearity  of  Phases. 
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as  a  Function  of  ©. 


Imaging  Interval  26  High  Non  -  Linearity 


03 
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Figure  64. 


Imaging  Interval  26  -  Poorly  Focused  with  High  Degree  of  Non- 
Linearity  of  Phases. 
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2.  Imaging  Intervals  with  2048  Pulses  and  4th-Degree  Motion  Compensa¬ 
tion 

Now,  if  4th-degree  motion  compensation  (i.e.,  motion  compensation  parameters 
{/21 5/22  5/23  j/24})  are  used,  the  following  non-linearity  of  phases  graph  is  generated  and 
shown  in  Figure  65: 

Degree  of  Non-Linearity  in  Each  Imaging  Interval 
6  -  Point  Delta  Wing  Data  Set  with  2048  Pulses  and  Degree  of  Motion  Compensation  -  4 


Figure  65.  Degree  of  Non-Linearity  of  Imaging  Intervals  of  2048  Pulses  and  4th- 

Degree  Motion  Compensation. 

In  general,  there  is  consistency  in  that  those  imaging  intervals  which  had  a  low 
degree  of  non-linearity  of  phases  in  Figure  58,  still  register  as  having  a  low  degree  of 
non-linearity  in  Figure  65.  The  same  applies  for  the  imaging  intervals  with  a  high  degree 
of  non-linearity.  There  are  two  reasons  for  the  differences  between  Figures  58  and  65. 
First,  the  AJTF  PSO  used  to  extract  the  search  parameters  has  a  certain  randomness  asso¬ 
ciated  with  its  search  pattern  and  it  is  possible  that,  for  certain  imaging  intervals,  incor¬ 
rect  search  parameters  were  extracted  and  thus  the  degree  of  non-linearity  calculation  was 
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slightly  off.  Secondly,  4th-degree  motion  compensation  is  more  accurate  so,  provided 
that  all  of  the  search  parameters  were  extracted  correctly,  it  provides  a  better  measure  on 
the  degree  of  non-linearity. 

Figure  65  indicates  that  the  best  three  imaging  intervals  are  imaging  intervals  15, 
14  and  23.  Imaging  intervals  15  and  23  are  the  same  as  before  and  are  shown  above. 
Imaging  interval  14  is  new  and  is  now  determined  to  be  one  of  the  best  three  imaging  in¬ 
tervals  in  the  data  set.  If  this  image  in  Figure  66  is  compared  to  imaging  interval  18  in 
Figure  61,  it  can  be  seen  that  imaging  interval  14  is  a  better  image  than  the  one  produced 
by  imaging  interval  18. 

Imaging  Interval  14  Low  Non  -  Linearity 


Figure  66.  Imaging  Interval  14  -  Well  Focused  with  Low  Degree  of  Non-Linearity 

of  Phases. 

Figure  65  indicates  that  the  three  worst  imaging  intervals  are  imaging  intervals 
26,  1 1  and  28.  Imaging  interval  12,  which  is  a  particularly  bad  imaging  interval  and  has 
been  shown  before  in  Figure  63,  falls  off  the  list  of  the  top  3  worst  imaging  intervals. 
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Fortunately,  it  is  still  ranked  as  having  a  high  degree  of  non-linearity  and  would  still  not 
be  imaged.  The  new  imaging  interval,  28,  is  shown  in  Figure  67  and  is  not  very  well  fo¬ 
cused  in  the  cross-range. 


Imaging  Interval  28  High  Non  -  Linearity 


Figure  67.  Imaging  Interval  28  -  Poorly  Focused  with  High  Degree  of  Non- 

Linearity  of  Phases. 

3.  Imaging  Intervals  with  1024  Pulses  and  2nd-Degree  Motion  Compen¬ 
sation 

There  are  several  reasons  to  use  imaging  intervals  with  fewer  pulses  in  the  cross¬ 
range.  The  most  obvious  is  that  is  that  the  TCP]  is  smaller  and  subsequently  there  should 

be  less  motion  error  in  each  imaging  interval.  Secondly,  since  there  are  more  imaging  in¬ 
tervals  available  for  processing,  if  an  imaging  interval  has  to  be  discarded,  there  are  still 
many  more  imaging  intervals  available  to  work  with.  The  downside  to  this  is  that  with 
less  pulses  in  the  cross-range,  the  cross-range  resolution,  A rcr ,  will  be  less. 
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Figure  68  shows  the  non-linearity  of  phase  graph  for  1024  pulses  in  the  cross¬ 
range  and  2nd-degree  motion  compensation.  The  three  best  imaging  intervals  are  imaging 
interval  45,  33  and  35  and  are  shown  below  in  Figures  69-71.  The  three  worst  imaging 
intervals  are  imaging  intervals  56,  25  and  29  and  are  shown  in  Figures  72  -  74.  Again, 
imaging  intervals  with  a  low  degree  of  non-linearity  are  easily  focused  while  imaging  in¬ 
tervals  with  high  degree  of  non-linearity  do  not  focus. 

Degree  of  Non-Linearity  in  Each  Imaging  Interval 
6  -  Point  Delta  Wing  Data  Set  with  1024  Pulses  and  Degree  of  Motion  Compensation  -  2 


Figure  68.  Degree  of  Non-Linearity  of  Imaging  Intervals  of  1024  Pulses  and  2nd- 

Degree  Motion  Compensation. 
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®ideal  ®ideal 


Figure  69.  Imaging  Interval  45  -  Well  Focused  with  Low  Degree  of  Non-Linearity 

of  Phases. 


Imaging  Interval  33  Low  Non  -  Linearity 
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Figure  70.  Imaging  Interval  33  -  Well  Focused  with  Low  Degree  of  Non-Linearity 

of  Phases. 
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Imaging  Interval  35  Low  Non  -  Linearity 
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Figure  7 1 .  Imaging  Interval  35  -  Well  Focused  with  Low  Degree  of  Non-Linearity 

of  Phases. 


Imaging  Interval  56  High  Non  -  Linearity 


®ideal  ®ideal 


Figure  72.  Imaging  Interval  56  -  Poorly  Focused  with  High  Degree  of  Non- 

Linearity  of  Phases. 
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Imaging  Interval  25  High  Non  -  Linearity 


®ideal 


Figure  73.  Imaging  Interval  25  -  Poorly  Focused  with  High  Degree  of  Non- 

Linearity  of  Phases. 


Imaging  Interval  29  High  Non  -  Linearity 
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Figure  74.  Imaging  Interval  29  -  Poorly  Focused  with  High  Degree  of  Non- 

Linearity  of  Phases. 
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4.  Imaging  Intervals  with  1024  Pulses  and  4th-Degree  Motion  Compensa¬ 
tion 

If  3D  motion  detection  is  repeated  for  1024  pulse  imaging  intervals  but  this  time 
using  4th-degree  motion  compensation,  the  degree  on  non-linearity  graph  is  generated  as 
in  Figure  75.  It  can  be  seen  from  this  graph  that  imaging  intervals  10,  33  and  2  possess  a 
low  degree  of  non-linearity  while  imaging  intervals  24,  56  and  29  possess  a  high  degree 
of  non-linearity.  Good  imaging  intervals  10  and  2  are  shown  below  in  Figures  76  and  77 
and  poor  imaging  interval  24  is  shown  in  Figure  78.  Again,  there  is  consistency  between 
the  2nd  and  4th-degree  linearity  of  phase  graphs  (in  Figures  68  and  75).  Imaging  intervals 
that  register  as  having  a  high  degree  of  non-linearity  of  phases  for  the  2nd-degree  motion 
compensation  also  register  as  having  a  high  degree  of  non-linearity  of  phases  for  4th- 
degree  motion  compensation. 

Degree  of  Non-Linearity  in  Each  Imaging  Interval 
6  -  Point  Delta  Wing  Data  Set  with  1024  Pulses  and  Degree  of  Motion  Compensation  -  4 


Figure  75.  Degree  of  Non-Linearity  of  Imaging  Intervals  of  1024  Pulses  and  4th- 

Degree  Motion  Compensation. 
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Imaging  Interval  1 0  Low  Non  -  Linearity 


Figure  76.  Imaging  Interval  10  -  Well  Focused  with  Low  Degree  of  Non-Linearity 

of  Phases. 


© 
4 — 

O 

£= 

o 

r-J  -43 

0  tt 


cn 

LU 

CD 


2 

1.5 

1 

0.5 


100  200  300 

© 

ideal 


— 

— 

— 

____  V 

- - 

- 1 

/ 

r - 

_ _ 

...V! 

_ 

L _ _ 

* 

...  2! 

V 

f , 

100  200  300 

© 

ideal 


_ v 

jf 

7 

/ 

_ j 

_ _ j 

/ 

/ 

✓ 

7 

^  1000 
© 

'S  800 


■tr-  o 
rn  -jzi 

©  t3 


600 

400 


w  200 
LU 

TO 


.7. 

y 

r 

... 

72.. 

. 

/ 

y 

7 

j* 

r 

100  200  300 

© 

ideal 


Figure  77.  Imaging  Interval  2  -  Well  Focused  with  Low  Degree  of  Non-Linearity  of 

Phases. 
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Imaging  Interval  24  High  Non  -  Linearity 
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Figure  78.  Imaging  Interval  24  -  Poorly  Focused  with  High  Degree  of  Non- 

Linearity  of  Phases. 


5.  Comparison  between  the  Different  Imaging  Interval  Sizes  and  the  De¬ 
gree  of  Motion  Compensation 

Now  that  the  degree  of  linearity  of  phases  has  been  generated  for  the  two  different 
imaging  interval  sizes  and  2nd  and  4th-degree  motion  compensation,  it  is  possible  to  com¬ 
pare  the  results  and  infer  which  combination  is  better  suited  for  rapid  ISAR  image  focus¬ 
ing.  The  most  straight-forward  way  to  go  about  this  is  in  the  form  of  a  table  that  lists  the 
processing  time  required  to  produce  each  graph.  Note  that  the  total  processing  time  is  the 
time  that  it  takes  to  bring  every  imaging  interval  in  the  60,000  pulse  6  -  Point  Delta  Wing 
experimental  data  set  from  their  non-focused  state  to  the  point  were  their  degree  of  3D 
motion  has  been  measured  and  a  decision  can  be  made  with  respect  to  whether  or  not  the 
imaging  interval  should  be  focused  or  discarded.  The  results  are  shown  in  Table  2: 
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#  of  Pulses  in 

Cross-range 

Number  of 

Images 

Degree 

Total  Process¬ 
ing  Time 

(seconds) 

Processing 

Time  per  Image 

(seconds) 

2048 

29 

2 

561.8 

19.3 

2048 

29 

4 

1710.1 

58.9 

1024 

58 

2 

547.8 

9.4 

1024 

58 

4 

1664.1 

28.7 

Table  2.  Processing  Time  -  Degree  of  Non-Linearity  of  Phases  for  Different  Im¬ 
aging  Intervals. 


The  first  thing  that  is  noticeable  from  Table  2  is  that,  unsurprisingly,  the  process¬ 
ing  time  for  the  entire  data  set  increases  by  over  a  factor  of  3  when  4th-degree  motion 
compensation  is  used  instead  of  2nd-degree  motion  compensation.  Since  there  is  only  a 
marginal  increase  in  the  quality  of  the  images  (only  imaging  interval  13  with  2048  pulses 
is  significantly  affected),  using  higher  order  motion  compensation  for  the  6  -  Point  Delta 
Wing  experimental  data  set  is  a  poor  use  of  processing  time. 

This  now  leaves  the  2048  pulse  imaging  interval  with  2nd-degree  motion  compen¬ 
sation  to  be  compared  its  1024  pulse  counterpart.  Using  the  data  set  divided  into  1024 
pulse  imaging  intervals  has  two  advantages  over  the  2048  pulse  imaging  intervals.  First, 
it  has  twice  as  many  images  to  work  with  so,  when  the  AJTF  PSO  does  not  converge  to  a 
focused  image  (which  will  occasionally  be  the  case  when  the  PSO  is  set  at  an  optimal 
combination  of  probability  of  convergence  and  run  time),  discarding  an  imaging  interval 
is  no  great  loss.  Secondly,  half  of  the  processing  time  per  image  is  helpful  in  getting  im¬ 
ages  to  the  ISAR  operator  as  quickly  as  possible.  The  data  set  with  2048  pulses  in  the 
cross-range  has  the  advantage  of  producing  better  images  with  point  scatterers  that  are 
well  separated  in  Doppler.  In  an  operational  situation,  these  trade-offs  will  have  be  bal¬ 
anced  against  each  other.  However,  for  the  purposes  of  this  thesis  it  was  easier  to  set  up 
3D  motion  detection  for  the  2048  imaging  intervals  since  3D  motion  detection  and  the 
AJTF  algorithm  are  very  sensitive  not  only  to  range  bin  selection  but  as  well  to  point 
scatterer  selection. 
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C.  CHAPTER  SUMMARY 

This  Chapter  presented  Thayaparan’s  method  of  3D  motion  detection  after  it  had 
been  adapted  to  the  AJTF  PSO.  It  was  very  effective  and  efficient  at  sorting  between  the 
good  imaging  interval  which  focused  well  and  the  poor  imaging  intervals  that  possessed 
3D  motion  and  could  not  be  focused.  The  poor  imaging  intervals  can  simply  be  dis¬ 
carded  and  not  presented  to  the  ISAR  operator.  In  addition,  it  was  shown  that  for  the  6  - 
Point  Delta  Wing  experimental  data  set,  2nd-degree  motion  compensation  is  the  most  effi¬ 
cient  used  of  computational  time. 
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VI.  CONCLUSION 


A.  SUMMARY  OF  RESULTS 

1.  AJTF  GA  and  PSO  Algorithms 

The  AJTF  GA  and  PSO  algorithms  performed  very  well  in  their  task  of  reducing 
the  computational  burden  required  by  the  AJTF  algorithm  when  employing  the  exhaus¬ 
tive  search  for  search  parameter  extraction.  Typically,  less  than  50%  of  the  cost  function 
calculations  were  required  for  the  GA  and  PSO  algorithms  to  extract  the  translational  and 
rotational  motion  compensation  parameters  which  were  able  to  focus  the  ISAR  images  in 
the  cross-range.  In  addition,  a  new  method  of  determining  search  parameter  suitability, 
the  automatic  recognition  of  focus  using  the  fast  Fourier  transform,  was  designed  and 
proved  very  fast  and  extremely  effective  in  providing  accurate  cost  function  returns  that 
gave  an  accurate  figure  of  merit  on  the  suitability  of  solution  space  search  parameters. 
When  comparing  the  two  different  algorithms,  the  PSO  algorithm  consistently  outper¬ 
formed  the  GA  and  was  able  to  converge  to  a  focused  image  with  less  calculations  of  the 
cost  function.  Also  of  important  note  is  that  the  use  of  the  PSO  algorithm  in  ISAR  image 
focusing  is  a  unique  application  of  this  evolutionary  search  technique  which  has  not  been 
published  in  the  scientific  literature. 

2.  3D  Motion  Detection 

Thayaparan’s  3D  motion  detection  algorithm  in  [12],  after  adaptation  to  this  the¬ 
sis’  AJTF  PSO  algorithm  and  uniquely  changing  the  measure  of  non-linearity  of  phases 
to  the  percentage  of  deviation  from  the  line  of  least-squares  fit,  proved  effective  in  sepa¬ 
rating  good  imaging  intervals,  possessing  only  2D  motion,  from  the  poor  imaging  inter¬ 
vals  that  possessed  a  significant  amount  of  3D  motion.  Application  of  the  AJTF  PSO  to 
the  3D  detection  algorithm  significantly  reduced  computational  time,  especially  since  it 
was  required  that  3  phase  functions  be  extracted  from  3  different  point  scatterers.  The 
AJTF  PSO  algorithm  is  a  unique  approach  to  the  3D  motion  detection  problem  and  cre¬ 
ates  an  algorithm  which  is  both  fast  and  accurate.  It  was  also  realized  that  the  applica¬ 
tion  of  motion  compensation  greater  than  the  2nd-degree  just  adds  to  the  computational 
burden  associated  with  3D  motion  detection  without  adding  any  significant  improvement 
to  3D  motion  detection  or  image  quality. 
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B.  FOLLOW-ON  WORK 

For  those  interested  in  ISAR  image  formation,  there  are  several  areas  of  research 
that  could  be  studied.  First,  an  important,  but  classified,  thesis  topic  could  involve  the 
AJTF  PSO  algorithm  being  tested  against  ISAR  returns  of  actual  military  targets  to  verify 
that  the  motion  compensation  algorithm  can  be  deployed  against  such  targets  of  interest. 
Secondly,  jet  engine  modulation  (JEM)  introduces  severe  time-variance  into  an  ISAR  re¬ 
turn.  Worthwhile  research  would  look  at  ways  of  removing  the  image  smearing  due  to 
JEM  or  classifying  targets  based  on  their  JEM  signature.  A  good  reference  for  this  re¬ 
search  topic  is  [15]. 
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APPENDIX 


MATLAB  CODE 


Listed  in  this  appendix  is  the  Matlab  code  used  in  this  thesis.  The  order  of  listing 
is  the  GA  code,  followed  by  the  PSO  code  and  the  3D  motion  detection  code.  The  code 
requires  MatLab  6.5  or  later  and  the  signal  processing  toolbox  to  run. 

A.  GA  CODE 


function  [fDl,  fD2,  FS,  FST]  =  GA_Image_Focus(G,  Rl,  R2,  m_orderl,m_order2,num_pop,num_gen, output) 
%  Written  by  Capt  Wade  Brinkman,  Canadian  Forces 
%  Oct  9,  2004 
%  Modified:  10  Nov  2004 

%  As  partial  fulfillment  of  the  MSEE  Degree  at  NPS 
%  Original  DRDC  version  written  by:  Jonathan  Waisman 
%  This  version  used  the  exhaustive  search  technique  as  described  in 
%  thesis  -  see  reference  page  for  DRDC  paper 


%  This  algorithm  will  focus  IS  AR  images  using  a  Genetic  Algorithm 
%  application  of  the  AJTF  algorithm. 


%  Input: 

%  G  -  signal  (2D  MxN  array) 

%  Rl  -  a  range  bin  that  will  be  used  for  translational  motion  comp  (  1  to  M  ) 
%  R2  -  a  range  bin  that  will  be  used  for  rotational  motion  comp  (1  to  M); 

%  m  orderl  -  order  of  motion  in  Rl 
%  m_order2  -  order  of  motion  in  R2. 

%  num_pop  -  the  number  in  the  GA  population 
%  num  gen  -  number  of  generations  in  the  run 
%  output  -  1  -  output  TF  Transforms 
%  0  -  supresses  the  TF  Transforms 

%  Output: 

%  fDl  -  the  doppler  coefficients  for  the  trans  motion  comp 
%  fD2  -  the  doppler  coefficients  for  the  rotational  motion  comp 
%  FS  -  final  focused  signal 
%  FST  -  signal  after  translational  motion  comp 


format  short; 

[M  ,  N]  =  size(G);  %M  -range  bins  N  -  pulses 


fDl  =  abs(fft(G(Rl,:))); 
fDl  =  find(fDl  ==  max(fDl)); 


fD2  =  abs(fft(G(R2,:))); 

ID  2  =  find(fD2  ==  max(fD2)); 


%find  the  factorials  for  l..m_orderl  &  two 
for  k  =  1  :m_orderl ; 

Orderlfac(k)  =  factorial(k); 
end 


for  k  =  1  :m_order2; 

Order2_fac(k)  =  factorial(k); 
end 


%look  at  the  uncompensated  signal; 
if  max(size(G))  ==256; 

srdi(G,Rl  ,R2, 1 ,256, 1 ,64, 1 );  colormap('hof  ); 
else 
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srdi(G,Rl  ,R2,round(N/2)-20,round(N/2)+20, 17,31 ,2); 
end 

%now  we  can  look  at  the  TF  of  the  two  range  bins 
if  output  ==1; 
figure; 

tfr_xl  =  tfrsp(G(Rl,:)');  %stft 

tfr_xl  =  [tfr_xl(N/2+l:N,:);  tfr_xl(l:N/2,:)]; 

imagesc(abs(tfr_x  1 )); 

title(['Uncompensated  -  Range  Bin  #  ',num2str(Rl)]); 
colormap('hof); 


figure; 

tfr_x2  =  tfrsp(G(R2,:)');  %stft 

tfr_x2  =  [tfr_x2(N/2+l  :N,:);  tfr_x2(l:N/2,:)]; 

imagesc(abs(tfr_x2)); 

title(['Uncompensated  -  Range  Bin  #  ',num2str(R2)]); 
colormap('hof); 
end; 


%%%%%%%%%%%%%%%%%%%%%0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o%%% 

%'Step  1 :  Translational  motion  compensation  using  the  first  range  bin' 


[fDl  FST]  =  GA_AJTF(G,Rl,m_orderl,num_pop,num_gen,0,Orderl_fac,fDl); 


%%%%%%%%%%%%%%%%%%%%%0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o%%% 

%%%%%% 

%'Step  2  -  Display  Translational  Motion  Compensation  Results' 
if  max(size(G))  ==  256; 

srdi(FST,Rl  ,R2, 1 ,256, 1 ,64, 1 );  colormap('hot'); 
else 

srdi(FST,Rl  ,R2,round(N/2)-20,round(N/2)+20, 17,31,1); 
end 


%display  the  2  range  bins  after  translational  rotation 
if  output  ==  1 

XI  =  FST(R1,:);  %lst  range  bin 

tfr_xl  =  tfrsp(Xl');  %stft 

tfr_xl  =  [tfr_x  1  (N/2+ 1  :N, : ) ;  tfr_x  1(1  :N/2,:)]; 

figure;  imagesc(abs(tfr_xl)); 

title(['Range  Bin  ',num2str(Rl),'  after  Translational  Motion  Compensation']); 
colormap('hot'); 


X2  =  FST(R2,:);  %lst  range  bin 
tfr_x2  =  tfrsp(X2');  %stft 
tfr_x2  =[tfr_x2(N/2+l  :N,:);tfr_x2(l  :N/2,:)]; 
figure;  imagesc(abs(tfr_x2)); 

title(['Range  Bin  ',num2str(R2),'  after  Translational  Motion  Compensation']); 
colormap('hot'); 
end; 


%'Step  3  -  Rotational  motion  compensation  using  the  second  range  bin' 


[fD2  FS]  =  GA_AJTF(FST,R2,m_order2,num_pop,num_gen,l,Order2_fac,fD2); 


%'Step  4  -  Display  Rotational  Motion  Compensation  Results' 
if  output  ==  1; 

XI  =  FS(R1,:);  %lst  range  bin 

tfr_x  1  =  tfrsp (XT);  %stft 

tfr_xl  =  [tfr_x  1  (N/2+ 1  :N, : ) ;  tfr_x  1(1  :N/2,:)]; 

figure;  imagesc(abs(tfr_xl)); 

title(['Range  Bin  ',num2str(Rl),'  after  Rotational  Motion  Compensation']); 
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colormap('hot'); 


X2  =  FS(R2,:);  %lst  range  bin 
tfr_x2  =  tfrsp(X2');  %stft 
tfr_x2  =[tff_x2(N/2+l  :N,:);tfr_x2(l  :N/2,:)]; 
figure;  imagesc(abs(tfr_x2)); 

title(['Range  Bin  ',num2str(R2),'  after  Rotational  Motion  Compensation']); 
colormap('hot'); 
end; 


%finally  the  compensated  image: 
if  max(size(G))  ==256; 

srdi(FS,Rl  ,R2, 1 ,256, 1 ,64, 1 );  colormap('hof ); 
else 

srdi(FS,Rl,R2,round(N/2)-20,round(N/2)+20, 1 7,3 1 ,2); 
end 


function  [fD,NS]  =  GA_AJTF(SIG,R,m_order,num_pop,num_gen,ToR,Order_fac,fD); 


%Capt  Brinkman,  Canadian  Forces 
%Oct  23,  2004 
%Modified:  20  Nov  04 

%this  function  is  main  Genetic  Algorithm  function  that  calls  all  the 
%other  GA  functions 

%The  goal  of  this  real  valued  GA  algorithm  is  to  converge  very  quickly  to 
%the  doppler  parameters  of  the  basis  function 


%SIG:  the  signal  to  be  analyzed; 

%R:  range  bin 
%m_order:  order  of  motion 

%num_pop:  this  is  the  population  of  our  genetic  pool 
%num_gen:  the  number  of  max  #  of  generations  to  run 
%ToR:  0  -  compensating  translational  motion 
%  1  -  compensating  rotational  motion 

%Order_fac  -  precalculated  factorials  needed  for  the  basis  function 
%fD  -  fl  search  parameter 


%Output  -  fD  -  returns  the  doppler  coefficients  for  the  motion  comp 
%  NS  -  motion  compensated  signal 


[rbins  pulses]  =  size(SIG); 
index  = 1 ; 

%Step  1 :  Find  the  initial  Population  and  rank  it 
population  =  I  POP(m  order,  num_pop, pulses, fD); 

%creation  our  initial  population  of  size  num_pop 
%with  values  from  -pulses :pulses 
[gen_fit]  =  FITNESS_FFT(population,SIG(R,:),ToR,Order_fac); 
%find  the  fitness  of  the  first  population 
mut_gen  =  num_gen; 
while  num_gen>0 


genscore(index)  =  sum(genfit); 

%ranks  the  overall  fitness  of  the  population 
generation_value(index)  =  sum(gen_fit)./num_pop; 

%sum  of  the  parents  raw  score  -  this  should  increase  over  time 
bestvalue(index)  =  max(genfit); 

%this  is  the  first  best  value  found; 


%Step  2:  Based  on  Suitability,  Generate  Children 
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children  =  CHILDS(population,gen_fit,gen_score(index)); 


%Step  3:  Mutation 

mut_prob  =  (0.3);  %30%  mutation  rate. 

children  =  MUTATION(children,mut_prob, pulses, morder); 


%Step  4:  Determine  Suitability  of  Children: 

[fit_child]  =  FITNESS_FFT(children,SIG(R,:),ToR,Order_fac); 


%Step  5:  Reinsertion 

%inserts  the  best  of  old  population  and  new  population  into  the  new 
%population 

[population, genfit]  =  REINSERTION (population,  children,  gen_fit,fit_child); 


numgen  =  numgen  -1; 
index  =  index  +1; 


end; 


%pick  out  the  best  function; 

best  =  population(fmd(gen_fit  ==  max(gen_fit)),:); 

fD  =  best(l,:); 

%Now  the  new  signal  can  be  calculated 
if  ToR  ==  0;  %translational  motion 
hp  =  basis2(fD,Order_fac,m_order,l, pulses, ToR); 
if  pulses<2%tum  off  for  now 
tfrhp  =  abs(stft(hp')); 

tfr_hp  =  [tfr_hp(pulses/2+l:pulses,:);  tfr_hp(l:pulses/2,:)]; 
figure;imagesc(abs(tfr_hp)) ; 
end; 

NS  =  SIG.*conj(repmat(hp,rbins,l)); 
end; 


if  ToR  ==  1;  %rotational  motion  comp 
clear  theta; 

[dummy,  theta]=basis2(fD,Order_fac,m_order,  1 , pulses, ToR); 
spacing  =  ((max(theta)-min(theta))/pulses)*0.999; 
uniform_samples  =  min(theta)+(l:pulses)*spacing; 

NS  =  interplq(theta.',SIG.',uniform_samples.').';  %reformat  the  radar  data 
%so  that  is  uniformly  samples  in  theta 

end; 


figure; 

plot(generation_value); 
hold  on;  plot(best_value,'g.-'); 


function  [population]  =  I_POP(m_order,num_pop,N,fD); 


%function  I  POP  will  create  the  initial  population  of  possible  solutions 

%m_order  -  order  of  motion 

%num_pop  -  max  number  in  the  population 

%N  -  the  values  of  the  pop  will  be  between  -N:N; 

%fD  -  is  our  fl  coefficient  as  determined  by  the  FFT 

%if  morder  ==  2; 

%  population  =  linspace(-N,N,num_pop)'; 

%elseif  m_order>2; 

%  population  =  [linspace(-N,N,num_pop)'  (-N  +  2*N*rand(num_pop,m_order-2))]; 
%end; 


population  =  -N  +  2*N*rand(num_pop,m_order-l); 
%guess  at  f2,  f3,  f4 . m  order  as  random  vars 
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fD_p(l:num_pop,:)  =  fD; 
population  =  [fD_p  population]; 
%add  fl  to  the  population 


function  [genfit]  =  FITNESS  FFT(pop,  X,  ToR,fac); 
format  long 

%Capt  Brinkman,  16  Oct  2004 
%Modified  20  Jan  2005 
%this  function  generates  the  suitability  of  each 
%parent  in  the  current  generation 


%pop  -  matrix  of  the  current  population 

%  -  fl,  f2,  f3  etc  depending  on  the  degree; 

%X  -  the  range  bin  under  examination 

%ToR  -  0  -  translational 

%  -  1  -  rotational 

%fac  -  precalculated  factorials 


%what  we  are  doing  here  is  trying  to  drive  the  slope 
%of  the  scatter  in  the  T  -  F  plane  to  zero.  When  this 
%  is  achieved,  translational  or  rotational  motion  can  be 
%removed. 


[num_pop,  degree]  =  size(pop); 


[a  pulses]  =  size(X);  %number  of  pulses  in  the  range  bin 


if  ToR  ==  0;  %compensating  translational  motion: 
hpvec  =  basis2(pop,fac, degree, num_pop, pulses, ToR); 

%hp_vec  is  the  basis  func  for  the  pop 
X  =  repmat(X,num_pop,l); 
test_result  =  X  .*  conj(hp_vec); 

test  =  (abs(fft(test_resulf)).A2)';%find  the  power  spectral  density  of  pop 
for  k  =  1  :num_pop; 
test_fft  =  test(k,:); 


max_test  =  max(test_fft); 

avg_PSD  =  max_test/sum(test_fft);  %Percentage  of  Power  under  max  point 


gen  fit  1  (k,  1  )=avg_PSD;  %max_sum; 

%when  the  percentage  of  energy  in  the  dominant  scatterer  of  the 
%range  bin  is  a  maximum,  we  have  our  best  focus 
end; 


gen_fit=  100.A(gen_fitl); 
gen_value  =  gen_fitl; 


theta  =  []; 
end; 


if  ToR  ==1 ;  %compensting  rotational  motion 

[dummy  theta]  =  basis2(pop,fac,length(fac),num_pop, pulses, ToR); 
for  k  =  1  :num_pop; 

spacing  =  ((max(theta(k,:))-min(theta(k,:)))/pulses)*0.999; 
uniformsamples  =  min(theta(k,:))+(l:pulses)*  spacing; 
test_sig  =  interplq(theta(k,:)',X',uniform_samples')'; 

%reformat  the  signal  so  that  it  is  linear  sample  w/ 

%theta  not  time. 

test_fft  =  abs(fft(test_sig)).A2;  %fmd  the  power  spectral  density  of  the  scatterer  rb 
max_test  =  max(test_fft); 

avgPSD  =  max_test/sum(test_fft);  %Percentage  of  Power  under  max  point 
gen_fit  1  (k,  1  )=avg_PSD ;  %max_sum; 

%when  the  percentage  of  energy  in  the  dominant  scatterer  of  the 
%range  bin  is  a  maximum,  we  have  our  best  focus 
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end;  %end  for  loop 
end; 


gen_fit=  100.A(gen_fitl); 

%the  exponential  provides  separation  between  good  values 
%and  bad  values.  The  higher  the  exponential,  the  greater  the 
%GA  searches  around  first  global  max  found.  The  smaller,  the 
%longer  it  takes  the  pop  to  base  towards  global  max  and  a 
%global  search  is  executed, 
genvalue  =  genfitl; 


function  [children]  =  CHILDS(population,gen_fit,gen_score); 
%Capt  Wade  Brinkman 
%Canadian  Forces  -  20  Oct  04 


%this  function  will  select  the  best  parents  based  on  their  assessed  fitness 
%and  breed  them  in  order  to  create  the  children; 


%a  roulette  wheel  selection  method  is  used.  This  is  discussed  in  detail  in 
%Goldberg  and  the  GA  toolbox  manual 


%population  -  the  current  GA  pop  that  we  are  working  with 
%gen_fit  -  score  between  1  and  zero  that  assesses  their  fitness 
%gen_score  -  the  score  for  the  current  gen  (sum(gen  fit)); 

[N  degree]=  size(population); 

%Step  1 :  Assign  Percentage  of  the  roulette  wheel 


rawscore  =  floor((gen_fit./gen_score).*  10000)71 0000; 

%converts  gen_fit  to  a  percentage  truncated  at  4  decimals; 
wheel  =  cumsum(raw_score); 

%the  values  progressively  increase  and  this  becomes  our  random 
%roulette  wheel 


%Step  2:  find  the  parents  chosen  by  the  roulette  wheel 
select  =  max(wheel)*rand(l,length(gen_fit)); 

%these  are  our  pointers  for  the  entire  selection  process 
select  =  repmat(select,N,l); 

%convert  to  matrix 

wheel  =  repmat(wheel,  1  ,length(gen_fit)); 

%convert  to  matrix  for  vector 


choice  =  wheel>=select;  %places  a  1  where  this  is  true 

%the  index  of  the  first  value  in  each  column  corresponds  to 
%a  parent  that  has  been  chosen  by  the  random  roulette  wheel 

for  index  =  1:N; 

check  =  find(choice(:,index)==l); 

Pool(index,:)=population(check(l),:);  %necessary  for  version  6 
%check  =  find(choice(:,index)==l ,  1 , 'first');  -works  for  version  7 
%Pool(index,:)=population(find(choice(:,index)==l,l, 'first'),:); 
%this  is  our  list  of  good  parents  that  will  produce  good  offspring 
end; 


pooll  =Pool(l:N/2,:); 

pool2  =  Pool(N/2+l  :N,:);  %break  the  pool  up  into  two  equal  matrix 


%the  first  one  in  the  population,  after  first  gen  best  parent,  is 
%guaranteed  two  spots  -  one  in  each  pool 
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pooll(l,:)  =  population^,:); 
pool2(l,:)  =  population(l,:); 


indexl  =  randperm(N/2); 
index2=  randperm(N/2); 
%generate  a  random  order; 


alpha  =[ones(N/2,l)  (0  +  l.*rand(N/2, degree- 1))];  %alpha  is  our  scaling  factor 
%alpha  =  repmat(alpha,l, degree);  %reformat  it  into  a  N/2  -  degree 
%array  to  allow  for  vectorization 
childl  =  alpha.*pooll  (indexl, :)+(l -alpha). *pool2(index2,:); 
child2  =  (1-alpha).  *pooll(indexl,:)+(alpha).*pool2(index2,:); 

%calculate  our  children  using  the  method  described  in 
%thesis  and  Junfei  Li  and  Hao  Ling's  paper  -  see  references 


children  =  [child  1  ;child2] ; 

%combine  our  children  array  and  return. 


function  [population, gen_fit,gen_value]  =  REINSERTION (population, children, genfit,... 

fit_child,gen_value,value_child); 

%Capt  Wade  Brinkman 
%Modified:  20  Nov  04 
%Input 

%population  -  old  population 

%children  -  the  children  of  the  current  population 

%gen_fit  -  fitness  of  each  element  in  population 

%fit_child  -  fitness  of  the  children 

%Output 

%population  -  new  population 

%gen_fit  -  new  generation  fitness 

%gen_value  -  new  generation  value 

%This  file  will  create  a  new  population  based  on  fitness 


%first  step  is  to  combine  the  two  vectors 
temp_pop  =  [population;  children]; 
tempfit  =  [gen _fit;fit_child] ; 


for  k=l  :length(population)/2; 

max  fit  =  max(temp_fit);  %find  the  maximum  fitness 
max_index  =  find(temp_fit  ==  max_fit);  %maximum  index 
max  index  =  max  index(l);  %ensures  one  value  is  found 


population(k,:)  =  temp  j)op(max_index, :); 
gen_fit(k,:)  =  temp_fit(max_index); 


temp_fit(max_index(l))  =  0; 
end 


index_new  =  find(temp_flt~=0); 

new_fit(  1 :  length(index_new))=temp  fit(index  new); 

new_pop(l:length(index_new),:)  =  temppop(index_new,:); 


for  k  =  (length(population)/2+l):length(population); 
index_rand  =  round(l+(length(population)-l)*rand); 
population(k,:)  =  new_pop(index_rand,:); 
gen_fit(k,:)  =  new_fit(index_rand); 
end; 


%this  is  our  new  population  with  our  new  measures  of  fitness; 


function  [show  image]  =  srdi(sig,rbl,rb2,xmin,xmax,ymin,ymax,itype); 
%si  =  srdi(FS,28,19,505,520,17,31,2);  will  work  for  fsmall  data  set 
rbinl  =  sig(rbl,:); 
rbin2  =  sig(rb2,:); 
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showimage  =  rdi(sig); 


subplot(2,2,[l  3]); 
imagesc(show_image); 
ylabel('Down  Range  -  Range  Bins'); 
xlabel('Cross  Range  -  Doppler  Bins'); 


axis([xmin  xmax  ymin  ymax]); 
if  itype  ==  0; 

title('Uncompensated  Image'); 
elseif  itype  ==  1 ; 

title('Image  After  Translational  Motion  Compensation'); 
elseif  itype  ==  2; 

title('Final  Image  after  Rotational  Motion  Compensation'); 
end 


frbinl  =  fliplr(abs(fftshift(fft(rbinl))).A2); 
f_rbin2  =  fliplr(abs(fftshift(fft(rbin2))).A2); 
n  =  2:(length(rbinl)+l); 


subplot(2,2,2); 

plot(n,f_rbinl,'r.-'); 

title(['Power  Spectrum  Range  Bin  ',num2str(rbl)]); 
axis([xmin  xmax  min(f_rbinl)  max(f_rbinl)]); 
ylabel(' Amplitude'); 
xlabel('Cross  Range  -  Doppler  Bins'); 


subplot(2,2,4); 

plot(n,f_rbin2,'r.-'); 

title(['Power  SpectrumRange  Bin  ',num2str(rb2)]); 
axis([xmin  xmax  min(f_rbin2)  max(f_rbin2)]); 
ylabel(' Amplitude'); 
xlabel('Cross  Range  -  Doppler  Bins'); 


function  [hp, theta]  =  basis2(pop,fac, degree, num_pop, pulses, ToR) 

%Capt  Wade  Brinkman 
%16  Oct  04 
%Modified  20  Jan  05 

%this  function  will  calculate  the  basis  function  for  our  current  population 

%the  calculation  will  be  done  in  parallel  to  speed  up  comp  time 

%this  version  removes  the  for  loops  which  cuts  down  the  run  time  by  half 

%this  is  important  since  we  are  trying  to  drive  the  run-time  of  the  GA  and  Swarm  to 

%almost  zero 

%Input: 

%pop  -  current  population 

%fac  -  the  factorials  used  in  the  calculation 

%degree  -  the  degree  of  motion 

%num_pop  -  the  number  in  current  population 

%pulses  -  #  of  pulse  in  cross  range 

%ToR  -  Translational  or  rotational  0  -  theta  =  [];  1  -  hp  =  []; 

%Output: 

%hp  -  the  basis  function  of  each  mbr  of  pop 
%theta  -  the  phases  associated  with  each  mbr  of  pop 
hp  =  zeros(num_pop,  pulses); 


t  =  ((l:pulses)./pulses)';  %t  =[(1  2  3  4  ...  pulses)./pulses]'  -  column  vector 
ttemp  =  t; 
for  i  =  2:  degree 
t  =  [t  (t_temp).Ai]; 

%tAl  forms  a  column  of  t,  tA2  forms  the  next... used  to  vectorize 
end 


pop  fac  =  pop./(repmat(fac,num_pop,l));  %predividing  the  population  by  its  factorials 
theta  =  (pop_fac(l:num_pop,:)*t(l:pulses,:)'); 


106 


if  ToR  ==  1; 

hp  =  []; 

else 

hp(l  :num_pop,l  :pulses)  =  exp(theta.*(-j*2*pi)); 
theta  =  []; 
end; 

%finally,  this  calculates  the  basis  function  for  the  entire  population, 
function  f  =  rdi(s) 

%  generates  a  range  doppler  image  of  the  signal  s 


f  =  fft(s'); 

[n,m]  =  size(f); 

f  =  [f(round(n/2)+l:n,:);  f(l:round(n/2),:)];  %  equivalent  to  fftshift 

f  =  abs(f/l); 
f  =(  f./max(max(f))); 

figure;  imagesc(abs(f)); 
colormap('jef); 


B.  PSO  CODE 

function  [FS  fDT  fDR]  =  Swarm_Image_Focus_3D(G,  R1,D  T,  RFocus, D  R, num_particles, cycles) 
%  Written  by  Capt  Wade  Brinkman,  Canadian  Forces 
%  Nov  13,  2004 
%  Modified:  12  Mar  2005 

%  As  partial  fulfillment  of  the  MSEE  Degree  at  NPS 
%  Original  DRDC  version  written  by:  Jonathan  Waisman 
%  The  DRDC  version  used  the  exhaustive  search  technique  as  described  in 
%  thesis  -  see  reference  page  for  DRDC  paper 
%  This  algorithm  will  focuse  IS AR  images  using  a  Swarm  Algorithm 
%  application  of  the  AJTF  algorithm. 

%  This  algorthm  uses  a  Swarm  Algorithm  as  described  in  Jacob  Robinson 
%  and  Daniel  Boeringer  papers  -  see  reference  page  -  with  modifications 
%  to  fit  the  AJTF  algorithm 

%  This  algorithm  uses  several  scatterers  from  the  same  range  bin  for 
%  rotational  motion  compensation. 

%  Input: 

%  G  -  signal  (2D  MxN  array) 

%  R1  -  a  range  bin  that  will  be  used  for  translational  motion  comp 
%  D_T  -  degree  of  translational  motion  compensation 
%  RFocus  -  [Range  Bin,  #  Scatterers,  Range  Bin,  #  Scatterers,....] 

%  list  of  range  bin  and  number  of  scatterers  for  rotational  motion 
%  compensation 

%  D  R  -  degree  of  rotation  motion  compensation 
%  num_particles  -  the  number  of  particles  in  the  Swarm  -  identical  to 
%  population 

%  cycles  -  number  of  times  the  Swarm's  position  are  updated  -  identical 
%  to  number  of  generations  in  GA 

%  Output: 

%  fDT  -  the  doppler  coefficients  for  the  translational  motion 
%  compensation 

%  fDR  -  the  doppler  coefficients  for  the  rotational  motion  compensation 
%  FS  -  final  focused  signal 


format  short; 

[M  ,  N]  =  size(G);  %M  -range  bins  N  -  pulses 
NRot  =  length(RFocus); 
if  mod(N_Rot,2)~=0; 

disp ('Rotational  Range  Bin  Matrix  incorrectly  formatted'); 
FS  =  G;fDT  =  [];  fDR  =  []; 
return 
end 
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N_RSB=  N_Rot/2;  %number  of  range  bins  for  rot  mot  comp 
RFocus  =  reshape(RFocus,2,N_RSB); 


for  k  =  1  :D_T;%find  the  factorials 
Orderlfac(k)  =  factorial(k); 
end 

for  k  =  1  :D_R; 

Order2_fac(k)  =  factorial(k); 
end 


%%%%%%%%%%%%%%%%%%%%%0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o%%% 

%'Step  1 :  Translational  motion  compensation  using  the  first  range  bin' 

%For  this  step,  the  Swarm_AJTF_3D  will  return  the  Doppler  Frequencies,  but  not 
%focus  the  image 


[fDT]  =  Swarm_AJTF_3D(G,Rl,D_T,num_particles,  cycles,  0,Orderl_fac); 

FST  =  tr  focus  3D(G,fDT,0); 

%%%%%%%%%%%%%%%%%%%0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o0/o%%% 


FST  t  =  FST; 

%'Step  2  -  Rotational  motion  compensation  using  the  Range  Bins 
%and  number  of  scatterers  found  in  RFocus 
r  fDR  =  1 ;  %index  for  the  rotational  doppler 
for  k  =  1  :N_RSB 

for  g  =  1  :RFocus(2,k); 

[fDR(r_fDR,  :)]=Swarm_AJTF_3D(FST_t,RFocus(  1  ,k), DR, num_particles, cycles,  1  ,Order2_fac); 
if  RF  ocus(2  ,k)>  1  &&g<RF  ocus(2  ,k) 

X_t  =  fft(FST_t(RFocus(l,k),:)); 
index  =  fDR(r  fDR,l)+(-3:3); 
for  kk  =  1  :length(index); 

if  index(kk)<l  X_t(N+index(kk))  =0; 

elseif  index(kk)>N  X_t(index(kk)-N)=0;  %remove  used  scatterer  from  Range  Bin 
else  X_t(index(kk))  =0;  end; 
end; 

FST_t(RFocus(l,k),:)  =  ifft(X_t); 
end; 

rfDR  =  rfDR  +1 ; 
end 
end; 


FD  =  mean(fDR,l); 

FS  =  tr  focus  3D(FST,FD,1); 

function  [fD]  =  Swarm_AJTF_3D(SIG, R,m_order,num_particles, cycles, ToR,Order_fac,fD); 


%Capt  Brinkman,  Canadian  Forces 
%Nov  13,  2004 
%Modified  20  Jan  05 

%this  function  is  main  Swarm  Algorithm  function  that  calls  all  the 
%other  Swarm  functions.  It  is  called  once  for  translational  and  once  for 
%rotational  motion  compensation. 

%The  goal  of  this  real  valued  Swarm  algorithm  is  to  coverge  very  quickly  to 

%the  doppler  parameters  of  the  basis  function 

%this  is  a  line  search  -  it  will  find  fl,  fix  fl,  find  f2,  fix  f2  find  f3 

%etc.  Necessary  since  solution  space  grows  by  NAm_order 

%SIG:  the  signal  to  be  analysed; 

%R:  range  bin 
%m_order:  order  of  motion 

%num_particles:  this  is  the  number  of  partilces  in  our  Swarm 
%cycles:  max  #  of  generations  to  run 
%ToR:  0  -  compensating  translational  motion 
%  1  -  compensating  rotational  motion 

%Order_fac  -  precalculated  factorials  needed  for  the  basis  function 
%fD  =  fl 

%Output  -  fD  -  returns  the  doppler  coefficients  for  the  motion  comp 
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fD  =  abs(fft(SIG(R, :)));  fD  =  find(fD  ==  max(fD)); 
[rbins  pulses]  =  size(SIG); 
numcycles  =  cycles; 

GFit  Old  =  0; 

current_f  =  2;  %current  parameter  being  searched 
Swarm=zeros(num_particles,  1 ); 

Swarm(  1  :num_particles,  1  )=[fD] ; 


while  current  f<=m  order 


%Step  1 :  Find  the  initial  Swarm  and  evaluate  fitness 

[Swarm, velocity]  =  I_Swarm(Swarm,m_order,  num_particles, pulses, current_f); 

%creation  our  initial  swarm  and  velocity 
%with  values  from  -pulses  :pulses  and  intial  velocities  from 
%-pulses:pulses 
%eval  the  fitness 

[partfit]  =  FITNESS  FFT(Swarm,SIG(R, :),ToR, Order  fac(l :  curren  tf) ,  currentf) ; 

LBest  =  Swarm;  %these  are  the  local  'bests'  found  by  each  particle.  For  the  first  cycle,  by  default  it  is 
%equivalent  to  the  first  Swarm 
LFit  =  part  fit;  %LFit  is  the  fitness  of  the  local  best 
GFit  =  max(LFit);  %this  is  the  maximum  global  fitness  that  was  found 

GBest  =  LBest(fmd(LFit==max(LFit)),:);  %this  is  current  global  best . 

GFit  =  GFit(l);  GBest  =  GBest(l,:); 
index  = 1 ; 

partscore(index)  =  sum(LFit)./num_particles; 

%ranks  the  overall  fitness  of  the  population 
best_value(index)  =  GFit; 

%this  is  the  current  global  maximum; 


bestcount  =  cycles; 

o_Best  =  GFit;  %these  two  variables  are  used  to  cause  the  algorithm  to  exit  if  no  improvement  has  been  made  after  x  -  cycles 

%figure; 

while  cycles>0 

%hold  on;  plot(index,LBest(:,2),'g+');hold  on;plot(index,GBest(l,2),'b>');hold  on;  plot(index,Swarm(:,2),'c.'); 
cycles  =  cycles  -1; 
index  =  index  +1; 


%Step  2:  Update  velocities: 

[velocity]=v_update(velocity, LBest, GBest, Swarm, num_particles,current_f, pulses); 


%Step  3 :  Update  positions 

[Swarm, velocity]=p_update(Swarm, velocity, num_particles, pulses); 


%Step  4:  determine  fitness  of  the  Swarm 

[part_fit]  =  FITNESS_FFT(Swarm,SIG(R,:),ToR,Order_fac(l  :current_f),current_f); 


%Step  5:  determine  LBest  and  GBest 

[LBest,  LFit,  GBest,  GFit]  =  swap_best(LBest,  LFit,  GBest,  GFit,  Swarm, part_fit,current_f); 


part_score(index)  =  sum(LFit)./num_particles; 

%ranks  the  overall  fitness  of  the  population 
best_value(index)  =  GFit; 

%this  is  the  current  global  maximum; 


end; 


if  GFit  Old  >=  GFit 

GBest(current_f)=0;  %if  the  Fitness  of  this  degree  is  less  than  the  previous,  set  this  parameter 
%equal  to  zero. 

else 
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GFit_01d  =  GFit; 
end 


current_f  =  current_f+l; 
cycles  =  num  cycles; 
Swarm=repmat(GBest,numparticles,  1 ); 


end; 


fD  =  GBest; 


function  [Swarm,  velocity]  =  I_Swarm(Swarm,m_order,num_particles,N,current_f); 
%Capt  Wade  Brinkman,  Canadian  Forces 
%14  Nov  04 
%Modified  12  Jan  05 

%function  I  Swarm  will  ceate  the  initial  particles  of  possible  solutions 
%and  assign  velocities  random  in  both  mag  and  direction  to  each  particle 


%  Swarm  -  current  Swarm  -  will  be  expanded 
%m_order  -  order  of  motion 

%num_particles-  max  number  in  the  particles  in  the  Swarm 
%N  -  the  values  of  the  pop  will  be  between  -N:N; 

%fD  -  is  our  fl  coefficient  as  determined  by  the  FFT 


Swarm  =[Swarm  (-N  +  2*N*rand(num_particles,l))]; 

%guess  at  f2,  f3,  f4 . m  order  as  random  vars 

velocity  =  (2*N)*rand(num_particles,l); 

%magnitude  of  the  velocites 
v_dir  =  rand(num_particles,l)-0.5; 

%direction  of  the  particle  +ve  or  -ve 
index  1  =  find(v_dir>=0); 
index2  =  find(v_dir<0); 
v_dir(indexl)=l; 
v_dir(index2)=- 1 ; 


velocity  =  velocity.  *v_dir; 

velocity  =  [zeros(num_particles,current_f-l)  velocity]; 
%the  fl  velocity  will  always  be  fixed  at  zero. 


function  [velocity]  =  v_update(velocity, LBest, GBest, Swarm, num_particles,m_order, pulses); 
%written  by  Capt  Wade  Brinkman,  Canadian  Forces 
%Created  14  Nov  2004 

%This  function  will  update  the  velocities  based  on  the  particles 
%current  velocity  and  the  pull  of  the  the  Local  Best  and  the  GBest 
%Input: 

%velocity  -  current  velocity  vectors 

%LBest  -  matrix  of  local  best  discoveries 

%GBest  -  the  best  discovery 

%Swarm  -  current  search  location  of  the  particles 

%num_particles  -  number  of  particles  in  the  swarm 

%m_order  -  order  of  motion 

%pulses  -  number  of  pulses  -  also  corresponds  to  max  velocity 
%Output: 

%velocity  -  updated  velocity  vector 


K(  1  :num_particles,  1  :m_order)=0.729;  %constriction  factor 
rhol  =2.05; 
rho2  =  2.05; 

rvecl  =  rand(num_particles,m_order); 
r_vec2  =  rand(num_particles,m_order); 
Global_Best=repmat(GBest,num_particles,  1 ); 


no 


velocity  =  K.*(velocity  +  2.05.*r_vecl.*(LBest-Swarm)+2.05.*r_vec2.*(Global_Best-Swarm)); 
%this  sets  the  velocity  vector  for  the  particles 

%Now,  check  to  make  sure  that  no  particle  exceeded  the  max  velocity; 

[R  C  V]  =  fmd(velocity>2*pulses); 
if  length(V)>0; 
for  k=l  :length(R) 

velocity(R(k),C(k))  =  2*pulses; 
end 
end 

[R  C  V]  =  fmd(velocity<-2*pulses); 
if  length(R)>0 
for  k=l:length(R) 

velocity(R(k),C(k))  =  -2*pulses; 
end 
end 

function  [Swarm, velocity]=p_update(Swarm, velocity  ,num_particles, pulses); 

%Programmed  by  Capt  Wade  Brinkman,  Canadian  Forces 
%14  Nov  2004 

%based  on  the  new  velocities,  the  positions  will  be  updated.  Any  particle 
%outside  the  solution  space  of  +/-  pulses  will  have  its  velocity  in  that 
%dimension  zeroed 
%Input: 

%Swarm:  The  particle  Swarm 

%velocity:  Particles  Current  velocity 

%num_particles:  Number  of  particles  in  the  Swarm 

%pulses:  number  of  pulses  in  signal  which  also  defines  the  bounds 

%  of  the  solution  space 


Swarm  =  Swarm  +  velocity; 

[R  C  V]  =  find(Swarm>pulses); 


if  length(V)>0;  %+vewall 
for  k  =  1  :length(V); 

velocity(R(k),C(k))=-pulses*0.1*rand;  %rebounding  wall 
Swarm(R(k),C(k))=pulses; 
end 
end 


[R  C  V]  =  find(Swarm<-pulses); 


if  length(V)>0;  %-vewall 
for  k  =  1  :length(V); 

velocity(R(k),C(k))=pulses*0.1*rand;  %rebounding  wall 
Swarm(R(k),C(k))=-pulses; 
end 
end; 


function  [LBest,  LFit,  GBest,  GFit]  =  swap_best(LBest,  LFit,  GBest,  GFit,  Swarm, part_fit,m_order); 
%Written  by  Captain  Wade  Brinkman,  Canadian  Forces 
%14  Nov  2004 

%This  function  will  update  the  local  and  global  bests  as  required 
%Inputs: 

%LBest  -  matrix  of  local  best  parameters 
%LFit  -  a  column  vector  of  the  fitness  of  the  local  bests 
%GBest  -  is  the  parameters  of  the  global  best 
%GFit  -  fitness  of  the  global  best; 

%Swarm  -  current  parameters  of  the  particle  swarm 
%part_fit  -  current  Swarm's  fitness 
%m_order  -  order  of  motion 

%Output:  LBest,  LFit,  GBest,  GFit  -  updated  parameters 


MaxSwarmfit  =  max(part_fit);  %the  best  of  the  current  particles 
MaxSwarm  =  Swarm(find(part_fit  ==Max_Swarm_fit),:); 
if  Max_Swarm_fit>GFit; 


ill 


GFit  =  Max_Swarm_fit(l); 
GBest  =  Max_Swarm(l,:); 
end; 


vecl  =  LFit>=part_flt;  %is  current  LBest  better  than  current  location 
vec2  =  part_fit>LFit;  %is  current  location  better  than  LBest 


LFit  =  LFit.*vecl+part_fit.*vec2;  %replaces  the  better  current  positions  and  LFit 

LBest  =  LBest.  *repmat(vecl,l,m_order)+Swarm.*repmat(vec2,l,m_order);  %places  the  better  position  in  Local  Best 


function  [Focused_Signal]  =  tr_focus_3D(signal,fD,ToR); 
%Capt  Wade  Brinkman 
%15  Nov  2004 

%this  is  a  utility  function  tha  will  perform  translational 
%or  rotational  motion  compensation  on  the  signal 


[rbins  pulses]  =  size(signal); 
for  k  =  1  :length(fD); 

o_fac(k)  =  factorial(k); 
end 

[hp  theta]  =  basis2(fD,o_fac,length(fD),l, pulses, ToR); 


if  ToR  =  0 

Focused_Signal  =  signal.  *conj(repmat(hp, rbins,  1)); 
else 

spacing  =  (max(theta)-min(theta))/pulses*0. 99999; 
uniform_samples  =  min(theta)+(l:pulses)*spacing; 
FocusedSignal  =  interplq(theta.', signal. ',uniform_samples.').'; 
end; 


C.  3D  MOTION  DETECTION  CODE 

function  [NL]=MDetect_3D(namedata,num_pulses,load_data,save_data, output, Num_S,ordl,ord2,cmovie); 
%function  [NL]=MDetect_3D('data',2Al  1,1, 0,1, 2, 2, 2,0); 

%Capt  Wade  Brinkman 
%Created  6  Feb  2005 
%Modified  12  March  2005 

%This  is  a  program  that  will  detect  the  presense  of  3D  motion 

%It  uses  the  EAJTF  Algorithm  {as  described  by  Thayaparan's  unpublished  work} 

%with  the  Swarm  algorithm  as  its  search  engine 

%Linearity  of  Phases  is  the  Average  Linearity  of  Phases  Method 

%{as  described  by  Thayaparan's  unpublishedd  work} 

%INPUT 

%  namedata  -  the  name  of  the  data  set  {without  the  imaging  interval 
%  number} 

%  num_pulses  -  number  of  IS  AR  pulses  in  the  cross-range 
%  -  either 2A10or2All 

%  load  data  -  if  1  -  fD  and  theta  already  exisit 
%  save  data  -  save  the  fD  and  theta 
%  output  -  output  all  images 

%  Num_S  -  number  of  Scatterers  in  first  range  bin  {default  1  in  second} 

%  ordl  -  order  of  motion  comp  -  Translational 
%  ord2  -  order  of  motion  comp  -  Rotational 
%  cmovie  -  create  a  movie 
%OUTPUT 

%  NL  -  Degree  of  Non-Linearity 
fflag  =  1 ; 

numimages  =  floor(60000/num_pulses); 


%first  set  the  directory  of  choices 
if  (num_pulses  ==  2A1  l)&&(ordl==2); 
cd  LP2048  2D; 

directory  =  cd;  cd  ..;  %sets  our  active  directory  where  fD  and  theta  are  stored 
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rb_file  =  'rb_focus_3D';  %this  is  the  file  that  contains  the  range  bins  used 
tempload  =  cat(2, directory, '/',rb  file); 

load(temp  load);  %the  rb_focus_3D  contains  the  range  bins  for  3D  motion  detection 
rb_focus  =  rb_focus_3D;  clear  rb_focus_3D; 
elseif  (num_pulses  ==  2A1  l)&&(ordl==4); 
cd  LP20484D; 

directory  =  cd;  cd  %sets  our  active  directory  where  fD  and  theta  are  stored 

rbfile  =  'rb_focus_3D';  %this  is  the  file  that  contains  the  range  bins  used 
tempload  =  cat(2, directory, '/',rb_file); 

load(temp  load);  %the  rb_focus_3D  contains  the  range  bins  for  3D  motion  detection 
rb_focus  =  rb_focus_3D;  clear  rb_focus_3D; 
elseif  (num_pulses  ==  2A10)&&(ordl==2); 
cd  LP 1024  2D; 

directory  =  cd;  cd  %sets  our  active  directory  where  fD  and  theta  are  stored 

rb  file  =  'rb_focus2';  %this  is  the  file  that  contains  the  range  bins  used 
tempload  =  cat(2, directory, '/',rb_file); 

load(temp  load);  %the  rb_focus_3D  contains  the  range  bins  for  3D  motion  detection 
rbfocus  =  rb_focus2;  clear  rb_focus2; 
elseif  (num_pulses  ==  2A10)&&(ordl==4); 
cd  LP1024  4D; 

directory  =  cd;  cd  %sets  our  active  directory  where  fD  and  theta  are  stored 

rb_file  =  'rb_focus2';  %this  is  the  file  that  contains  the  range  bins  used 
tempload  =  cat(2, directory, '/',rb  file); 

load(temp  load);  %the  rb_focus_3D  contains  the  range  bins  for  3D  motion  detection 
rb  focus  =  rb_focus2;  clear  rb_focus2; 
else 

disp('Must  choose  2A1 1  or  2A10  pulses  or  Order  of  Motion  Error  equal  to  2  or  4'); 
return 
end; 


for  k  =1  :num  images; 

if  load  data  ==  0;  %need  to  generate  new  fDs  and  phases 

if  fflag  ==  1;  fflag  =  0;  fac  =  zeros(l,ord2);  for  g  =  l:ord2;  fac(g)  =  factorial(g);  end;  end;  %precalc  factorials 
image_int  =  cat(2,directory,'/',namedata,int2str(k-l));  %create  string  containing  name  of  imaging  interval 
load(image  int);  %load  imaging  interval  into  work  space 

[FS  fDT  fDR]  =  Swarm_Image_Focus_3D(dat,rb_focus(k,l),ordl,[rb_focus(k,2),Num_S,rb_focus(k,3),l],... 

ord2,20,100); 

if  output==  1 
rdi(FS); 

title(['Focused  Image  Imaging  Interval  ',num2str(k-l),'  -  ',num2str(num_pulses),'  Pulses  in  Cross  Range']); 
xlabel('Cross  Range  -  Doppler  Bins'); 
ylabel('Down  Range  -  Range  Bins'); 
end 

[dum  ord2]  =  size(fDR);  f  =  [fDR];  [m  n]  =  size(f); 

[hp, theta]  =  basis2(f,fac,ord2,m,num_pulses,l); 
elseif  load_data  ==  1 ;  %data  already  exisits 
doppler  =  cat(2,directory,'/','fD',int2str(k-l)); 
load(doppler); 

[r  c]  =  size(fD); 
fDT  =  fD(l ,:); 
fDR  =  fD(2:r,:); 
if  output  ==  1 ; 

imageint  =  cat(2,directory,'/',namedata,int2str(k-l)); 
load(imageint); 

FST  =  tr_focus_3D(dat,fDT,0); 

FS  =  tr_focus(FST,mean(fDR,  1),  1); 

title(['Focused  Image  Imaging  Interval  ',num2str(k-l),'  -  ',num2str(num_pulses),'  Pulses  in  Cross  Range']); 
xlabel('Cross  Range  -  Doppler  Bins'); 
ylabel('Down  Range  -  Range  Bins'); 
end; 

[dum  ord2]=size(fDR);  %detect  number  of  scatterers  and  order 
phases  =  cat(2, directory,'/', 'theta', int2str(k-l)); 
load(phases); 

if  fflag  ==  1;  fflag  =  0;  fac  =  zeros(l,ord2);  for  g  =  l:ord2;  fac(g)  =  factorial(g);  end;  end;  %precalc  factorial 
f  =  [fDR];  [m  n]  =  size(f); 
end; 

[NL(k)]  =  EDEGOFNONLINEARIT  Y2(  theta  ,f,min(size(theta)),fac,0  ); 
if  savedata  ==  1 ; 

doppler  =  cat(2,directory,'/','fD',int2str(k-l)); 
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fD  =  [fDT;fDR]; 
save(doppler,'fD'); 

phases  =  cat(2,  directory,  '/Vtheta',int2str(k-1)); 
save(phases,  'theta'); 
end; 
end; 


NL_norm  =  NL; 

figure;  plot(0:(num_images-l),NLjiorm,'r+-','linewidth',2);  grid  on;  axis  tight 
cutoff_line(l:(num_images))  =  mean(NL); 

%hold  on;  plot(0:(num_images-l),cutoff_line,'b-.','linewidth',2);  hold  off; 
xlabel('Imaging  Interval'); 
ylabel('Degree  of  Non-Linearity'); 

title(  {'Degree  of  Non-Linearity  in  Each  Imaging  Interval',... 

['6  -  Point  Delta  Wing  Data  Set  with  ',num2str(nurnj)ulses),... 

'  Pulses  and  Degree  of  Motion  Compensation  -  ',num2str(ordl)]}); 


%now  we  can  display  the  good  images 
%and  make  a  movie  if  the  option  is  selected 
[r  c  v]  =  fmd(NLnorm<=.  010*  mean  (NL)); 
for  k  =  l:length(c); 

imageint  =  cat(2,directory,'/',namedata,int2str(c(k)-l)); 
load(imageint); 

doppler  =  cat(2,directory,'/','fD',int2str(c(k)-l)); 
load(doppler); 

[row  col]  =  size(fD); 
fDT  =  fD(l,:); 
fDR  =  fD(2:row,:); 
fD  =  mean(fDR,l); 

FST  =  tr_focus_3D(dat,fDT,0); 

FS  =  tr_focus(FST,fD,l); 
image  =  rdi2(FS); 
bound  =  image>.3; 

[r2  c2  v2]  =  find(bound==l); 
centerr  =  round((max(r2)+min(r2))/2); 
center_c  =  round((max(c2)+min(c2))/2); 
axis([center_c-15  center_c+15  8  38]); 

if  cmovie  ~=1;  title(['Imaging  Interval:  ',num2str(c(k)-l)]);  end; 
if  cmovie  ==1; 
axis  off; 
pause(2); 

g_°4_NF(k)  =  getframe(gcf); 
end 
end 


if  cmovie  ==  1 ;  save  g_o4_NF  g _o4_NF;  end; 


%now  we  are  going  to  plot  out  the  three  best  and  three  worst  imaging 

%intervals 

NF_t  =  NF; 

for  k  =  1:3  %find  the  three  best 
temp  =  find(NF_t  ==  min(NF_t)); 
best_image(k)  =  temp-1 ;  NF  t(temp)  =  NaN; 
end 

for  k  =  1:3  %find  the  three  worst 
temp  =  find(NF_t  ==  max(NL_t)); 
worst_image(k)  =  temp-1;  NL_t(temp)  =  NaN; 
end 

%now  plot  out  the  3  best  and  3  worst  along  with  their  non-linearity  plots 


for  k  =  1:3 
figure; 

imageint  =  cat(2, directory, 7', namedata,int2str(best_image(k))); 
load(imageint); 
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phases  =  cat(2, directory,'/', 'theta', int2str(best_image(k))); 
load(phases); 

doppler  =  cat(2, directory, '/','fD',int2str(best_image(k))); 
load(doppler); 

[row  col]  =  size(fD); 
fDT  =  fD(l,:); 

fDR  =  fD(2:row,:);  fD  =  mean(fDR,l  ); 

FST  =  tr_focus_3D(dat,fDT,0); 

FS  =  tr_focus_3D(FST,fD,  1); 
imint  =  rdi2(FS); 
subplot(2,2,l); 
imagesc(imint); 

title(['Imaging  Interval  ',num2str(best_image(k)),'  Low  Non  -  Linearity']); 
bound  =  im_int>.3; 

[r2  c2  v2]  =  find(bound==l); 
center_r  =  round((max(r2)+min(r2))/2); 
center_c  =  round((max(c2)+min(c2))/2); 
axis([center_c-15  center_c+15  8  38]); 

[test]  =  EDEGOFNONLINEARITY 2(  theta  ,fDR,min(size(theta)),fac,  1  ); 
end; 


for  k  =  1:3 
figure; 

image_int  =  cat(2, directory, '/',namedata,int2str(worst_image(k))); 
load(imageint); 

phases  =  cat(2, directory,'/', 'theta', int2str(worst_image(k))); 
load(phases); 

doppler  =  cat(2, directory, '/VfD',int2str(worst_image(k))); 
load(doppler); 

[row  col]  =  size(fD); 
fDT  =  fD(l,:); 

fDR  =  fD(2:row,:);  fD  =  mean(fDR,l); 

FST  =  tr_focus_3D(dat,fDT,0); 

FS  =  tr_focus_3D(FST,fD,  1); 
imint  =  rdi2(FS); 
subplot(2,2,l); 
imagesc(imint); 

title(['Imaging  Interval  ',num2str(worst_image(k)),'  High  Non  -  Linearity']); 
bound  =  im_int>.3; 

[r2  c2  v2]  =  find(bound==l); 
center_r  =  round((max(r2)+min(r2))/2); 
center_c  =  round((max(c2)+min(c2))/2); 
axis([center_c-15  center_c+15  8  38]); 

[test]  =  E  DEG  OF  NONLINEARITY 2(  theta  ,fDR,min(size(theta)),fac,  1  ); 
end; 

function  [AVGN]  =  E_DEG_OF_NONLINEARITY2(  P  ,  fD  ,  L  ,fac,OP) 

%  AVGN  =  E  DEG  OF_N ONLINE ARITY 2  (  P  ,  fD  ,  L  )  - 
%  Finds  the  presence  of  3D  motion  in  each  imaging  interval.  And  returns 
%  the  following. 

% 

%  INPUT(S): 

%  P  -  the  phase  function  of  an  imaging  interval 
%  fD  -  doppler  coefficients  for  this  phase  function  (P  can  be  generated 
%  from  fD) 

%  L  -  number  of  point  scatterers  chosen  for  analysis 

% 

%  OUTPUT(S): 

%  AVGN  -  degree  of  3D  motion  for  the  imaging  interval  that  defined  the 
%  phase  function  P.  Note  that  AVGN  is  the  average  degree  of 
%  non-linearity  (or  in  other  words  the  degree  of  3D  motion)  for  the  non-linearities 
%  for  every  pair  of  phase  functions  defined  in  P. 

% 

%  Author  :  Jonathan  Waisman 

%  Modified  by  Capt  Wade  Brinkman,  Canadian  Forces 
%  First  Modified:  9  Feb  05 
%  Last  Modified:  29  Mar  05 


[H  ,  N]  =  size(  P  ); 

[B  ,  ORDER]  =  size(  fD  ); 
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Nij  =  zeros(  1  ,  L  ); 
AVGN  =  0; 

AProj  =  0; 


%ideal  doppler  coef  will  be  average 
ifD  =  mean(fD,l); 

%ideal  phase  function 

[hp  i_theta]=basis2(i_fD,fac, ORDER,  1  ,N,  1); 
aa  =  max(max([i_theta;P])); 
for  i  =  1  :L 

%  generate  the  non-linearity  between  the  IdealP  and  P(i,:) 

[Nij(  i  )]  =  EDEGOFNONLINEARITY 1  ( i_theta  ,  P(i,:),OP,i+l  ); 
%  sum  the  non-linearities  for  averaging 
AVGN  =  AVGN  +  Nij(  i ); 
end 


AVGN  =  AVGN  /  L; 

function  [  N12  Proj]  =  E  DEG  OF  NONLINEARITY  1  (i_P  ,  r_P, OP, index ) 

%  [  N12  ]  =  DEG  OF  NONLINEARITY  1  ( i  P  ,  r  P  )  -  calculates  the  degree  of 
%  non-linearity  between  i_P  and  r_P  using  a  least  squares  approach. 

% 

%  INPUT(S): 

%  r_P,i_P  -  2  phase  functions 

% 

%  OUTPUT(S): 

%  N12  -  degree  of  nonlinearity  between  these  2  particular  phase  functions 

% 

%  Author  :  Jonathan  Waisman 
%  Modified  by  Capt  Wade  Brinkman 
%  Last  Modified:  29  March  05 

%  i_P(t)  =  a  *r_P(t)  +  n(t)  is  the  relation  of  2  point  scatterers. 


p=polyfit(i_P,r_P,  1); 
p  =  fliplr(  p  ); 


N  =  length(  i  P  ); 
k  =  length(  p  ); 
func  =  zeros(  1  ,  N); 
t  =  linspace(l,max(i_P),N); 
for  x  =  1  :N 
for  i  =  1  :k 

func(  x  )  =  func(  x)  +  p(i)*(  i_P(x)  A  ( i  -  1  )  ); 
end 
end 


%  p  =  [p(l),p(2)]  where  least  squares  approximation  func(t)  =  p(l)  +  p(2)*t 
%  but  func  is  the  best  fit  linear  approximation  to  i_P,  and  thus: 

%  the  degree  of  nonlinearity  is  the  difference  in  the  plot  of 
%  plot(i_P,func)  and  plot(i_P,r_P); 


N12  =  0; 


P  =  r  P; 

FUNC  =  func; 
for  x  =  1  :N 

N12  =  N12  +  abs(P(x)-FUNC(x))./abs(FUNC(x)); 
end 


if  OP  ==  1 

subplot(2,2,index);  hold  on; 
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plot(i_P,func,'g.-Vlinewidth',2);hold  on; 
plot(i_P,r_P,'r-.71inewidth',2); 
xlabel('\Theta_i_d_e_a _f); 

ylabel({['\Theta_',num2str(index-l),'(t)'];'  as  a  Function  of  \Theta_i_d_e_a_l'}) 
grid  on;  axis  tight 
end; 


%to  see  the  degree  of  non-linearity,  enter  this  line  into  the  command  line: 
%figure;  plot(i_P,func,'g.-71inewidth',2);hold  on; 
%plot(i_P,r_P,'r.-71inewidth',2); 


%  The  degree  of  3D  motion  as  the  deviation  from  a  linear 
%  relationship  between  theta  and  phi  over  the  dwell  intervall 
%  is  directly  related  to  N  by  a  constant  Beta.  Thus  the  degree  of 
%  nonlinearity  will  give  a  degree  of  3D  motion! 
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